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CHAPTER  I 


INTRODUCTION 

During  the  Gulf  War,  Iraq  was  thought  to  possess  some 
30  fixed  SCUD  sites.  The  problem  with  these  sites  was  that 
Saddam  never  once  launched  a  SCUD  out  of  one  of  those.  He 
preferred  using  mobile  platforms  consisting  of  a  truck  that 
could  set  up,  launch  a  SCUD,  and  escape  in  a  matter  of 
minutes  [1] .  Even  though  they  were  militarily  insignificant, 
political  pressure  forced  many  forces  to  be  almost  solely 
utilized  for  search-and-destroy  missions  against  these 
targets.  This  presented  a  significant  problem  because 
western  Iraq  is  about  the  size  of  Vermont,  New  Hampshire, 
and  Massachusetts  put  together  [1]  . 

Needless  to  say,  a  better  job  could  have  been  done  in 
force  utilization  if  the  targets  could  have  been  found  more 
efficiently.  The  problem  then  becomes:  how  can  past 
movements  of  mobile  SCUD's  be  correctly  modeled  to  predict 
future  movements  and  narrow  search  areas  in  future 
conflicts?  This  will  decrease  the  time  and  amount  of 
resources  needed  to  destroy  SCUD's.  Fortunately,  that  is 
possible  in  the  future  if  space  assets  are  used  wisely. 

These  assets  can  provide  the  essential  initial  inputs 
of  environment  data  and  past  observations  needed  to  predict 
future  movements.  Satellites  like  Canada's  Radarsat  can 
provide  excellent  moisture  gradients  using  synthetic 
aperture  radar  (SAR)  [2] , [3] .  Moisture  content  is  important 


2 


because  mobile  Scud's  can't  operate  in  too  damp  of  soil. 
Other  space-based  assets  can  also  provide  remote  sensing 
capabilities  for  topographical  vector  fields  of  terrain. 
Aggregating  this  information  will  provide  a  basis  with  which 
to  work  on  predicting  future  movements  of  mobile  platforms 
like  SCUD  launchers. 

Predicting  future  movements  of  mobile  launchers  can 
hopefully  be  described  by  the  use  of  hidden  Markov  models 
(HMM) .  The  hidden  Markov  model  has  found  most  use  in  the 
area  of  speech  recognition  for  computers  but  has  recently 
found  use  in  several  different  areas.  Everything  from 
target  classification  to  computational  biology  now  uses  HMM 
[4] , [5] , [6] , [7] , [8] .  This  study  will  attempt  to  use  a  novel 
approach  to  HMM.  Based  on  artificial  observational  data 
generated  by  a  computer^,  an  algorithm  will  be  formulated  to 
determine  the  crew's  mindframe  at  any  time  iteration,  the 
correct  model  order,  and  make  predictions  for  possible 
future  movements.  The  objective  of  the  mobile  SCUD  is 
unknown,  so  the  prediction  of  future  movement  has  to  come 
from  past  observed  data.  A  full  definition  of  the  hidden 
Markov  model  can  be  seen  as  part  of  chapter  II. 

This  approach  is  novel  for  two  reasons.  The  first 
reason  has  to  do  with  the  time  of  analysis.  This  approach 
will  determine  useful  information  based  on  very  few 
observations.  There  is  no  prior  knowledge  that  the 
formulation  will  generate  any  useful  information  based  on  as 


1  No  actual  observations  of  bearing  were  available  to  use  as  data  for  this  study. 
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few  as  three  observations.  Also,  because  of  the  requirement 
for  timely  prediction,  it  is  important  that  the  initial 
parameters  of  the  model  be  accurate  based  on  apriori 
knowledge.  This  will  hopefully  reduce  the  observations  that 
are  required  to  get  information  about  the  object. 

The  second  novelty  deals  with  how  the  model  order  is 
chosen.  Past  observations  will  be  split  into  two  groups. 

The  data  designated  "distant"  past  will  be  used  to 
parameterize  the  model.  The  model  will  then  try  to  forecast 
the  data  designated  "recent"  past.  The  best  model  order 
will  forecast  the  recent  past  with  the  highest  probability 
[9] .  Verification  or  order  is  knowing  that  a  particular 
model  predicts  the  recent  past  with  a  higher  probability 
than  all  other  models.  In  other  words,  the  order  that 
predicts  the  recent  past  best  will  be  considered  the  proper 
order.  No  truth  model  is  designated  outright  in  this  paper 
to  generate  data.  Data  will  be  generated  from  an  ideal 
model  through  the  computer. 

It  is  important  to  note  that  available  data  will  be 
restricted.  With  such  small  amounts  of  data,  a  truth  model 
can  mislead  an  observer  over  a  small  amount  of  time  due  to 
its  probabilistic  nature.  Verification  of  model  order  under 
these  circumstances  would  be  difficult.  Therefore,  a  truth 
model  is  not  used  for  verification.  The  ideal  model  is  one 
that  would  predict  the  recent  past  with  perfect  accuracy. 

The  order  verification  will  come  in  how  close  the  given 
model  is  to  an  ideal  model,  not  the  truth  model. 
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Although  this  paper  emphasizes  mobile  SCUD  launchers  as 
an  application,  there  is  no  reason  why  this  algorithm  can't 
be  applied  to  other  objects.  For  that  reason,  instead  of 
saying  mobile  SCUD  throughout  the  paper,  we  will  refer  to 
what  is  being  observed  as  the  object. 

Scope  of  the  Study 

In  considering  how  to  predict  where  a  given  object  will 
be  at  some  point  in  the  future,  the  definition  of  the 
problem  is  critical  because  there  are  infinite  ways  to 
approach  the  problem.  The  design  has  to  reflect  the  desired 
outcome  and  initial  inputs.  The  desired  outcome  is 
information  made  up  of  bearing  sequence  predictions  and 
expected  protocols.  Movement  can  be  inferred  from  that 
information  depending  on  the  application.  The  observational 
data  will  come  in  the  form  of  direction  bearings.  Combined 
with  terrain  data,  these  initial  inputs  will  parameterize 
the  model . 

This  problem  will  begin  with  all  the  given  terrain 
information  assumed.  The  includes  road  maps,  weather, 
moisture  gradients,  vegetation,  etc.  This  is  significant 
because  a  lot  of  work  goes  into  getting  this  information. 

SAR  imagery,  road  maps,  and  terrain  maps  can  all  be  attained 
by  satellite.  The  scope  of  this  study  does  not  include  the 
positioning  of  these  satellites  at  the  proper  time,  remote 
sensing  any  of  the  images  to  a  specific  degree  of  accuracy, 
or  processing  these  images  into  usable  vector  fields.  This 
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study  will  start  with  the  terrain  information  and  proceed 
[3]  . 

An  important  assumption  is  that  all  the  prior 
activities  can  be  done  correctly  prior  to  this  analysis. 
Operationally  speaking,  there  are  three  different  kinds  of 
problems  that  can  occur:  modeling  error,  process  noise,  and 
measurement  noise.  The  first  error  happens  when  the  model 
order  is  incorrect.  Process  noise  comes  from  the  fact  that 
even  if  the  model  order  is  correct,  the  detail  is  not  the 
same  as  the  actual  system.  Lastly,  measurement  noise  comes 
from  the  instruments  or  the  environment  not  being  conducive 
for  accurate  measurement  [10] .  Assuming  that  the 
information  given  is  correct  nullifies  process  and 
measurement  noise.  Reducing  these  uncertainties  is  beyond 
the  scope  of  this  paper  so  we  assume  them  to  be  nonexistent. 
This  algorithm  will  try  to  select  the  model  order  and, 
therefore,  reduce  the  modeling  error.  That  will  lend  itself 
to  better  predictions  of  future  movements. 

After  analysis,  what  is  the  final  product?  Once  enough 
observations  of  the  object  are  taken,  the  expected  Markov 
state  at  any  time,  parameters  of  the  model,  and  predictions 
for  future  movements  will  be  the  output  information.  The 
Markov  state  of  the  object  will  correspond  to  what  state  of 
mind  the  object  is  in,  or  what  protocols  it  must  follow. 

The  parameters  of  the  model  will  include  how  the  object 
transitions  between  Markov  states  and  what  the  probability 
of  seeing  a  specific  movement  in  a  specific  state  is. 


6 


Finally,  the  prediction  of  future  movements  will  be  based  on 
the  parameters  of  the  model  and  will  show  trends  of  where 
the  object  is  most  likely  to  move  next. 

The  intent  is  to  present  an  algorithm  to  model  behavior 
of  an  object  with  HMM,  As  a  result,  no  actual  data  has  been 
used.  All  observations  are  a  consequence  of  a  random  number 
generator.  The  observational  data  used  for  validation  was 
designed  and  can  be  seen  as  part  of  Appendix  I .  There  is  no 
guarantee  that  the  object  will  behave  operationally  in  a  way 
that  lets  us  draw  meaningful  information  from  their 
movements  in  the  same  fashion  that  we  will  with  the 
generated  data,  but  that  is  an  unavoidable  uncertainty. 

In  order  for  this  model  to  approach  absolute  certainty, 
infinitely  many  observations  have  to  be  taken  over  time.  In 
an  operational  situation  this  might  not  always  be  possible. 
That  further  limits  the  ability  of  a  position  prediction  and 
order  determination  for  the  object.  In  other  words,  the 
extent  search  areas  may  be  narrowed  is  totally  dependent  on 
the  observation  time  allowed. 

Furthermore,  only  predictions  of  object  bearing  are 
given.  In  order  to  really  predict  future  movements  a 
velocity  and  time  interval  between  observations  would  have 
to  be  introduced.  This  is  dependent  on  the  characteristics 
of  the  object  and  the  type  of  orbit  used  for  observation, 
respectively. 

Trending  future  movements  is  not  an  exact  science. 
Highest  probable  sequences  will  be  given  based  on  observed 
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data  and  it  will  be  up  to  the  observer  to  draw  actual 
conclusions  about  bearing.  It  will  be  simple  in  most  cases 
but  that  is  entirely  dependent  on  the  nature  of  the  data  and 
how  much  time  the  object  is  observed. 

Understanding  the  mathematics  behind  this  formulation 
is  important  because  it  gives  insight  into  what  can  actually 
be  done  with  the  model .  The  methodology  of  the  algorithm  is 
the  most  important  part  of  this  study  because  of  the  new 
application  of  HMM.  Knowledge  of  the  methodology  is 
necessary  to  make  changes  to  the  algorithm  further  along  in 
the  paper.  Finally,  algorithm  validation  and  possible 
future  uses  of  the  formulation  will  be  presented. 


CHAPTER  II 


MATHEMATICS  BACKGROUND 

An  understanding  of  some  basic  math  concepts  behind  the 
hidden  Markov  model  is  essential  to  understanding  the 
definition  and  formulation  of  the  problem. 

Stochastic  Processes 

All  of  the  modeling  done  in  the  paper  can  be  understood 
within  the  confines  of  a  stochastic  process.  Generally 
speaking,  a  process  is  stochastic  if  for  every  X  event  in 
the  sample  event  space  S  a  real -valued  time  function  X(t,X) 
is  assigned.  For  each  event,  X,  there  is  a  single  time 
function  X  (t ,  A,j_)  =xi  (t)  called  the  sample  function.  The 
totality  of  sample  functions  for  a  system  is  called  an 
ensemble.  For  a  fixed  time  t  j ,  X(tj,A,)=Xj  is  a  random 
variable  [11] .  Figure  1  illustrates  the  notation.  Notice 
that  the  initial  conditions  of  the  sample  functions  are 
different  for  each  sample  function.  Also,  the  initial 
conditions  do  not  dictate  the  entire  sample  function. 

Stochastic  processes  can  also  be  defined  as  a 
collection  of  random  variables  {Xt-},  where  t  runs  through  a 
given  index  T.  Often  the  set  T  is  taken  in  time  so  it  is 
usually  non-negative  integers  [12] .  X^  is  some  measurable 
characteristic  at  each  t  in  the  index  T.  For  a  fixed  t  and 
X,  X (ti, =xi (ti)  will  become  a  number  [11] . 
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An  example  of  a  stochastic  process  is  the  toss  of  a 
coin  or  the  measurement  of  the  temperature  each  hour  of  the 
day.  In  the  case  of  the  coin  toss,  an  event  X,  might  be 
getting  three  heads  in  a  row.  The  sample  function  would 
look  like  three  heads  in  a  row.  The  random  variable  Xj  is 
the  side  of  the  coin  shown  at  the  particular  time  t 
consistent  with  the  sample  space.  Care  must  be  taken  when 
defining  stochastic  processes  because  probabilities  are 
assigned  to  events  and  not  outcomes.  In  the  case  of 
tracking  object  movements,  the  random  variable  can  be 
thought  of  as  the  measurement  of  a  object's  vector  position 
at  any  time  t  consistent  with  the  definition  of  the  sample 
space.  The  event  space  would  include  all  possible  bearing 
sequences  up  to  a  prescribed  time  t. 

A  numerical  outcome  whose  value  can  change  from 
experiment  to  experiment  can  be  thought  of  as  a  random 
variable.  The  formal  definition  is  a  real  valued  function 
whose  domain  is  a  sample  space  [13] .  From  the  temperature 
example,  the  sample  space  would  include  all  the  possible 
sequences  of  temperatures  for  that  season  up  to  a  time  t. 
The  real  valued  time  function  X{t,X)  would  be  the 
measurement  of  those  temperatures  sequences.  At  each  time 
step  the  random  variable  temperature  can  be  measured  from 
the  available  temperature  range  for  that  season. 

The  random  variables  are  assumed  discrete  for  this 
study.  The  random  variable  X  is  said  to  be  discrete  if  it 
can  take  only  a  finite  number  of  possible  values  of  X  [13] . 
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For  the  methodology,  an  outcome  of  the  model  is  assumed  to 
be  a  multiple  of  90  degrees.  More  complexity  requires  the 
outcome  multiple  to  be  smaller. 

Stochastic  processes  can  vary  in  the  form  they  take  on 
but  for  this  study,  at  particular  points  of  time  t  labeled 
0,  1,  ...,  the  model  is  found  in  one  Markov  state. ^  These 

states  are  mutually  exclusive  and  exhaustive  and  are  labeled 
0,  1,  ...,  M.^  The  points  in  time  are  spaced  depending  on 

the  behavior  of  the  physical  system  that  the  process  is 
embedded  in  [12] .  For  example,  a  measurement  can  be  taken 
at  a  change  in  the  system  or  at  regularly  spaced  intervals. 
For  our  purposes,  measuring  a  subject's  location  at  spaced 
intervals  is  most  appropriate  because  a  satellite  will  pass 
overhead  in  certain  time  intervals. 

Markov  Processes 

There  are  many  stochastic  processes  that  are  of 
interest.  A  Markov  process  is  just  a  stochastic  process 
that  exhibits  the  Markovian  property: 

=j%  =i}=JVi,^  =j\x,  =1}  , , , 


This  is  true  for  time  t  =  0,1,...  and  every  sequence  kQ, 
k^,...,  kt-i,  i,  j.  This  is  one  of  the  assumptions  that 
allow  stochastic  processes  analytical  tractability .  Other 

^There  are  only  a  finite  number  of  states  to  the  model. 

^Mutually  exclusive  means  that  the  elements  of  a  set  in  a  sample  space  do  not  overlap.  Exhaustive  means 
that  the  list  contains  all  possible  outcomes  [13]. 
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assumptions  about  the  joint  distribution  of  Xq,  X^, . . .  are 
also  necessary  and  will  discussed  later  in  the  paper  [12] 

Equation  (1)  basically  says  that  the  probability  of  a 
future  event  only  depends  on  the  present  state.  The 
conditional  probability  that  the  system  will  visit  any  other 
state  is  independent  of  all  states  except  the  present  state 
[12,13]  This  is  important  because  all  of  the  other  states 
that  the  model  has  visited  are  not  necessary  for  future 
predictions.  Essentially  this  assumption  makes  the  model 
simpler  to  define. 

In  the  model  dealing  with  location  prediction  the 
Markov  states  will  correspond  to  the  protocols  and  states  of 
mind  the  launch  team  has  to  operate  on.  In  other  words,  the 
Markov  state  will  correspond  to  the  rules  they  have  been 
given  and  what  frame  of  mind  they  are  in.  Different  states 
can  be  designed  to  mimic  different  possible  sets  of  orders. 
For  this  case,  the  Markovian  property  is  assumed.  The  crews 
operate  solely  on  their  current  orders  and  what  they 
perceive  danger  to  be  at  the  present  time.  Once  their 
Markov  state  changes  (i.e.,  their  orders  or  perception  of 
danger)  they  will  move  according  to  that  new  state 
exclusively. 


'*Let  Xj  and  X2  be  discrete  random  variables.  The  joint  probability  distribution  of  Xj  and  X2  is  given 
p(xi,X2)  =  P(Xi  =  xi,  X2  =  X2)  [13], 

A  and  B  are  any  two  events,  then  the  conditional  probability  of  A  given  B  is  given  by  P(A  |  B)  = 
P(AB)/P(B).  This  is  provided  that  P(B)>0.  Additionally,  if  two  events  are  independent  then  P(A  |  B)  = 
P(A)  [13], 
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The  notation  for  the  probability  of  going  from  one 
state  to  another  is : 


(«) 

y  (2) 

It  is  the  conditional  probability  that  the  random  variable 
X,  starting  at  state  i  will  be  in  state  j  after  exactly  n 
steps  [12] .  By  definition,  the  conditional  probabilities 
have  to  be  non-negative.  They  have  to  obey  the  following 
rules : 


>  0 

for  all  i  and  j;  n=0,l,2... 
M 

14”’=! 

for  all  i;  n=0,l,2,... 


(3) 

(4) 


These  rules  spring  from  the  fact  that  the  system  must  make  a 
transition  into  some  other  Markov  state  [12] . 

Ergodic  Systems 

The  ensemble  was  defined  above  as  all  the  sample 
functions  grouped  together.  It  is  often  useful  to  describe 
these  groups  of  sample  functions  in  the  context  of 
statistical  averages.  Consider  the  stochastic  process  X(t) . 
At  time  ti,  X(ti)=Xi  is  a  random  variable  with  distribution: 


Fx(^i'A)  = 


(5) 
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where  is  any  real  nuniber.  The  first  order  density 
function  is: 


The  mean  of  X(t)  can  be  found  by: 


(6) 


00 

MxUi)  =  E[X{t^)]  =  ^xfx{x,A)dx 

(7) 


If  the  stochastic  process  is  wide-sense  stationary,  the  mean 
is  constant  in  time  [11] . 


The  time-averaged  mean  of  a  sample  function  x(t)  of  a 
stochastic  process  X(t)  is  defined  as: 


1  ^ 

X  =  {x{t))  =  ^  x(,t)dt 


L 

1 


(9) 


By  being  stationary,  the  expected  value  of  both  sides  of 
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equation  (9)  can  be  taken.  It  will  yield  [11] : 

T_ 

E[x]  =  I  E[x{t)]dt  = 

2  (10) 

If  the  result  from  equation  (10)  is  assumed  to  be  true, 
the  time  averaged  mean  from  equation  (9)  and  the  ensemble 
mean  from  equation  (7)  are  equal.  By  definition  then,  an 
ergodic  process  is  one  that  has  equal  time  averages  for  each 
sample  function  and  are  also  equal  to  the  corresponding 
ensemble  averages  [11] .  The  key  to  a  process  being  ergodic 
is  to  be  wide- sense  stationary  and  have  statistics  that  do 
not  evolve  with  time  as  per  equation  (8) . 

This  result  may  seem  obscure  but  has  broad 
implications.  For  example,  if  a  process  is  ergodic,  then  by 
observing  a  single  sample  function  all  the  relevant 
statistics  can  be  determined  without  observing  more  sample 
functions  [11] .  In  other  words,  observing  the  process 
during  any  period  of  time  will  allow  a  full  description  of 
that  process  throughout  time.  This  process  is  assumed  to  be 
ergodic. 


Hidden  Harkov  model 

The  hidden  Markov  model  is  the  next  step  in  the 
complexity  of  the  Markov  process.  In  a  Markov  process,  the 
shift  from  one  Markov  state  to  another  is  probabilistic 
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while  the  output  is  deterministic.  The  hidden  Markov  model 
extends  this  simple  chain  to  a  doubly  stochastic  process. 

In  other  words,  both  the  output  symbol  and  Markov  state 
transition  are  determined  probabilistically.  The  sequence 
of  the  state  transitions  is  "hidden"  because  the  observer 
only  physically  sees  the  symbol  outputs  [4] ,  [14] ,  [15] ,  [16] . 

An  example  of  a  simple  Markov  process  is  a  coin  toss. 
The  process  is  given  two  states .  The  states  are  represented 
by  a  heads  or  tails  [14] .  Each  state  will  output  either  a 
heads  or  tails  as  the  observation  symbol,  respectively.  The 
output  symbol  has  no  probability  distribution  for  a  given 
state.  A  hidden  Markov  model  of  the  coin  toss,  on  the  other 
hand,  could  also  have  two  Markov  states.  These  states, 
however,  are  represented  by  differently  biased  coins.  Each 
coin  will  turn  up  either  heads  or  tails  based  on  a  different 
probability  distribution.  The  output  of  each  state  now  has 
a  probability  distribution  associated  with  it.  Each  state 
can  output  either  a  head  ot  tail .  The  process  has  become 
doubly  stochastic.  Modeling  coin  tosses  this  way  could  be 
important  because  the  more  complex  hidden  Markov  model  could 
be  better  suited  to  modeling  outputs. 

Because  the  hidden  Markov  model  (HMM)  is  a  doubly 
stochastic  process,  it  is  more  difficult  to  define.  Let  N 
be  the  number  of  Markov  states  in  the  model .  The  individual 
Markov  states  will  be  denoted  S={Si,  S2,...,  Sjj}  .  The 
number  of  distinct  observation  symbols  per  Markov  state  is 
given  by  M.  The  individual  symbols  are  denoted  V={vi, 
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V2,-../  ,  The  state  transition  probability  distribution 

A={aij}  is  given  by: 

;  1  <  i  and  j  <  N  (11) 

The  observation  symbol  probability  distribution  in  state  j , 
B={bj (k) }  is  given  by: 


Of  =  P[qM  =  Sj 


b,(k)  =  PMt)\q,=Sj] 


;  l<j^  and  l<k^  (12) 

Finally,  the  initial  state  distribution  n  =  {71^}  is  given 
by: 

;  l^i^  (13) 


With  the  parameters  of  N,M,A,B,  and  n  the  HMM  can  generate  a 
sequence  of  observations  0  =  Oi,  O2,...  O-p  [14]  ,  [16]  . 

Given  these  parameters,  several  problems  with  real-life 
applications  can  be  solved.  The  most  basic  problem  involves 
finding  an  efficient  way  to  compute  P(o|X),^  This  is  simply 
the  probability  of  observing  a  sequence  of  events  given  the 
model  parameters  [14]  .  The  highest  probability  sequence 
should  indicate  where  something  will  move  in  the  future. 
Examining  several  of  the  highest  probability  sequences, 
trends  should  be  able  to  be  drawn  as  to  where  objects  will 
be  later  in  time. 


denotes  the  parameters  of  the  model  A,  B,  and  n 
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The  second  problem  involves  finding  a  state  sequence,  Q 
=  ^1#  that  best  explains  the  observations  taken. 

This  commonly  involves  a  dynamic  programming  technique 
called  the  Viterbi  algorithm  [17] .  Although  no  totally 
correct  sequence  can  be  found,  it  does  show  most  probable 
sequences  and  can  lend  insight  into  the  model  [14] .  By 
knowing  what  state  the  teams  are  in,  we  can  lend  further 
insight  as  to  where  the  object  may  be  based  on  what  we  think 
their  protocols  are. 

Finally,  the  third  problem  involves  re-estimating  the 
parameters.  This  is  done  in  order  to  maximize  P(o|A,)  .  In 
the  end,  much  of  the  modeling  would  be  worthless  if  there 
was  no  way  to  improve  on  how  we  estimate  what  will  happen  in 
the  future  [14] . 

As  practical  as  this  type  of  modeling  seems,  it  does 
have  setbacks.  Numerically,  as  time  goes  to  infinity 
underflow  problems  for  the  probabilities  become  more  of  an 
issue.  Ironically,  a  infinite  number  of  samples  is  required 
to  generate  a  perfect  set  of  parameters  for  the  model 
[14] , [16] .  Secondly,  the  necessity  for  finite  training  sets 
forces  some  observations  to  0  probability  as  the  re¬ 
estimation  procedure  alters  the  observation  matrix 
[14] , [16] .  Finally,  the  re-estimation  procedure  tends 
toward  the  local  maximum  of  a  likelihood  function.  The 
likelihood  function  can  be  complex  with  many  local  maximums 
so,  estimating  initial  parameters  is  important  when  trying 
to  reach  a  global  maximum  and  speed  up  run  time  [14]  .  All 
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of  these  problems  will  be  addressed  later  in  the  methodology 
of  the  problem. 

Theoretically,  the  Markovian  property  itself  sometimes 
oversimplifies  the  nature  of  the  model.  Only  being 
dependent  on  the  current  state  is  sometimes  too  big  of  a 
simplifying  assumption.  Also,  HMM  is  only  well  defined  for 
single  independent  variable  function.  This  is  much  more  of 
a  problem  than  the  current  state  assumption  because  it  is 
currently  impossible  to  define  a  HMM  for  more  than  one 
independent  variable.  In  principle,  however,  HMM  can  be 
made  dependent  on  more  than  its  current  state  [4] . 


CHAPTER  III 


METHODOLOGY 

Referring  back  to  chapter  II,  the  HMM  needs  some  basic 
definitions  in  order  to  solve  the  problems  mentioned.  The  X 
parameters  need  to  be  defined,  with  the  number  of  Markov 
states,  and  the  number  of  observation  symbols  allowed. 

Also,  the  sequence  of  events  needed  to  reach  the  end  result 
has  to  be  presented. 

The  first  operation  after  estimating  the  initial 
parameters  and  gathering  observations  consists  of  finding 
the  forward  and  backward  probabilities.  These  probabilities 
help  to  find  probabilities  of  sequences  and  re-estimate 
parameters  of  the  model . 

Once  the  final  parameters  of  the  model  have  been  found 
through  re-estimation,  the  probability  of  the  most  likely 
sequences  can  be  found.  This  can  then  be  compared  to  what 
an  ideal  case  would  look  like.  The  particular  model  can  be 
evaluated  next  to  several  other  models  to  determine  the  best 
order.  Once  the  final  model  is  selected,  trends  can  be  seen 
as  to  where  the  object  will  go  and  what  its  protocols  are. 

Define  the  Model 

The  order  of  the  model  can  be  chosen  at  random.  The 
definition  of  the  model,  includes  the  number  of  states,  N, 
and  the  number  of  observations  allowed,  M.  The  methodology 
shall  start  with  a  fundamental  model  and  .work  to  more 
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complex  models  to  see  if  they  model  behavior  any  better.  To 
help  understand  the  significance  of  the  order  it  is  useful 
to  make  M  =  4  and  N  =  3,  initially.  Later  on,  there  will  be 
different  orders  to  consider  but  this  initial  model  is 
simplest  to  explain  and  understand  in  terms  of  observables 
and  Markov  states  visited.  The  4  observations  possible 
correspond  to  the  4  points  of  the  compass  (North,  East, 

South,  West) .  This  will  give  look  angles  of  90  degrees 
around  each  direction  of  the  compass.  In  the  MathCad 
formulation,  the  directions  of  the  compass  will  be 
represented  as  0,  1,  2,  and  3,  respectively.  The  3  states 
allowed  correspond  to  the  state  of  mind  of  the  team  driving 
the  launcher  and/or  the  rules  of  engagement  they  must 
follow. 

One  state  represents  trying  to  be  expedient.  The  team 
might  have  orders  to  get  out  of  the  area  as  soon  as  possible 
or  get  to  a  launch  site  by  a  certain  time.  Another  state 
represents  trying  to  be  concealed.  Orders  might  include 
hiding  in  the  local  area  or  evading  incoming  warplanes .  The 
last  state  represents  one  of  ambivalence.  Lacking  all 
instruction  and  information,  traveling  in  any  direction  will 
be  equally  likely. 

These  definitions  are,  of  course,  not  the  only  way  to 
define  models.  The  number  of  observations  we  allow  is  only 
limited  by  computer  speed.  The  included  look  area  that  each 
observation  has  will  only  decrease  with  more  observations 
allowed.  The  number  of  states  allowed  is  limited  by 
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computer  speed  and  creativity.  Other  than  what  was 
mentioned  above,  states  can  be  also  defined  by  moving  versus 
stationary  orders,  pre-launch  versus  post-launch  status,  or 
any  combination  thereof.  Through  our  testing  with  MathCad, 
the  number  of  states  was  the  most  significant  factor  in 
determining  run  time  on  the  computer.  The  scope  of  this 
study  will  focus  on  4  and  6  observations  with  3,  4,  and  5 
states  allowed.  Adding  states  by  re- formulating  the  model 
will  be  an  exercise  to  examine  the  order  of  the  model  and 
will  be  estimated  with  a  beta  density  function. 

Estimate  initial  matrices 

In  order  to  begin  the  simulation  an  initial  model  has 
to  be  defined.  These  matrices  represent  apriori 
probabilities  of  the  model.  With  these  probabilities,  the 
model  can  be  started  close  to  what  the  optimal  model  is. 

The  better  the  algorithm  is  in  estimating  the  matrices,  the 
faster  the  re-estimation  will  be. 

The  most  difficult  matrix  to  estimate  was  the 
observation  probability  matrix,  B.  The  MathCad  algorithm 
can  be  seen  as  part  of  Appendix  A.  With  all  terrain, 
weather,  and  road  maps  at  hand,  the  individual  vector  fields 
can  be  aggregated  to  form  one  final  field.  As  mentioned 
above,  this  operation  is  in  the  realm  of  image  processing 
and  assumed  given.  Therefore,  we  can  estimate  an  aggregate 
field  right  away  and  proceed  with  the  algorithm.  A  possible 
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example  of  this  aggregate  can  be  seen  as  part  of  Figure  2 . 

It  is  a  relative  plot  to  show  what  the  vector  field  might 
look  like.  Units  and  components  will  be  specified  by  the 
application.  It  is  assumed  that  the  higher  sections  of  the 
z-function  correspond  to  poor  weather  conditions,  bad  roads, 
or  good  concealment.  Conversely,  launch  minded  expedient 
teams  will  tend  to  move  towards  the  lower  sections  of  the  z- 
function.  This  convention,  however,  will  also  be  chosen  to 
fit  the  application. 

Another  convention  chosen  was  to  begin  the  object  in 
the  middle  of  the  aggregate  field.  Then  by  differencing 
that  element  with  every  other,  the  object  could  start  at  a 
relative  ground  zero.  The  values  of  the  new  vector  field 
can  be  seen  as  part  of  matrix  M  on  the  page  following  the 
contour  plot  in  Appendix  A. 

The  differencing  above  poses  a  problem.  The  M  matrix 
can  have  negative  elements  if  the  middle  element  is  higher 
valued  than  other  elements .  Negative  numbers  are  not 
acceptable  because  these  will  translate  later  into 
probabilities  that  would  be  infeasible.  Here,  this  problem 
was  remedied  by  forcing  a  small  positive  number  on  any 
negative  sum.  Of  course,  problems  can  arise  from  such  a 
simple  approach  (i.e.,  starting  at  the  top  of  a  hill)  but, 
it  will  suit  fine  for  most  applications. 

Now,  by  looking  in  each  direction  the  maximum  value  can 
be  determined  by  looking  at  the  associated  look  angle. 

Adding  the  maximum  z-function  values  together  in  a  weighted 
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Figure  2.  3-D  Plot  of  Aggregate  Fields 
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sum  across  the  look  area  will  translate  into  the  apriori 
probabilities  for  that  direction. 

These  algorithms  can  be  seen  as  part  of  Appendix  A. 

They  are  named  according  to  the  direction  in  which  they 
look.  The  bars  associated  with  each  direction  denote 
MathCad's  programming  structures.  The  east  program,  for 
example,  finds  the  maximum  z- function  value  of  the  matrix  M 
looking  in  the  east  direction.  Care  must  be  taken  in 
deciding  where  to  look  because  the  top  of  the  matrix  doesn't 
correspond  to  north  on  the  map  from  figure  2.  The  0,0 
element  on  figure  2  is  at  the  lower  left  corner  while  the 
0,0  element  on  the  matrix  is  in  the  upper  left. 

For  the  east  direction,  look  to  the  top  of  the  M 
matrix.  Starting  at  the  right  and  working  left  in  columns'^, 
each  applicable  row  is  examined  and  the  highest  value  is 
stored®.  A  while  loop  is  given  to  each  column  in  the  look 
area  for  that  direction.  A  weighting  is  given  to  each 
column  according  to  how  directly  it  faces  the  corresponding 
direction.  The  column  directly  east  is  weighted  the  highest 
while  columns  on  the  fringe  of  the  look  area  are  given  a 
smaller  weight.  The  west  direction  is  the  same  as  east  but 
will  look  down  the  matrix  M  instead.  With  the  other  two 
directions,  the  rows  will  take  the  place  of  columns  and  vice 
versa.  As  the  look  areas  become  smaller  due  to  more 
allowable  observations,  the  number  of  while  loops  will 

^The  columns  are  denoted  the  local  variable  b  in  these  programming  structures 

®A  row  is  considered  applic^le  if  it  falls  within  that  direction.  A  row  towards  the  bottom  of  the  matrix 

would  not  have  anything  to  do  with  being  in  the  east  direction. 
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decrease  while  the  number  of  total  programming  structures 
will  have  to  increase. 

Dividing  each  direction  by  the  sum  of  all  directions 
will  normalize  the  sum  to  unity.  Negative  numbers  in  this 
operation  would  cause  meaningless  probabilities.  Depending 
on  how  the  states  are  defined  for  a  particular  model,  these 
numbers  can  be  manipulated  to  represent  the  probable 
directions  for  that  state.  The  straight  forward  calculation 
above  makes  for  acceptable  estimates  of  the  likelihood  of 
moving  in  a  particular  direction  for  a  conceal  state,  for 
example .  This  is  true  because  of  the  way  the  convention  was 
arranged  earlier.  If,  on  the  other  hand,  the  object  needed 
to  move  quickly,  staying  on  the  low  portions  of  the  vector 
field  would  be  desirable.  This  can  be  described  by  taking 
the  reciprocals  of  the  direction  numbers  and  adding  them 
together  the  same  way  to  get  unity.  For  a  complete 
description  of  how  each  Markov  state  calculated  its  initial 
observational  symbol  likelihood's  see  chapter  IV  dealing 
with  these  descriptions . 

Whatever  the  formulation,  the  rows  of  the  B  matrix  have 
to  equal  the  number  of  states  allowed  and  the  summation  of 
each  row  has  to  be  unity  as  explained  above.  The  number  of 
columns  of  the  B  matrix  has  to  equal  the  number  of 
observation  symbols  allowed. 

The  state  transition  matrix.  A,  is  not  part  of  the 
physical  data.  This  matrix  shows  how  the  object  in  question 
transitions  between  operating  modes  (i.e.,  Markov  states). 
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The  formulation  will  initially  assume  that  the  object  will 
most  likely  transition  back  to  the  state  it  left  from  in  one 
step.  In  other  words,  the  changes  in  state  typically  won't 
occur  because  the  Markov  state  is  relatively  stable  over 
time . 

The  beta  density  function  will  be  used  to  estimate  the 
transitions  from  state  to  state  of  the  HMM.  The  beta 
density  function  is  convenient  to  use  because  it  is  easy  to 
parameterize  to  allow  for  recurrence  and  it  terminates  at 
unity.  Segmenting  the  density  function  according  to  the 
number  of  states  and  integrating  will  give  starting  points 
for  the  A  matrix  that  make  physical  sense  in  accordance  with 
the  assumption  above. 

Appendix  A  experiments  with  several  different  methods 
for  representing  the  matrices.  They  serve  as  a  check  on  the 
most  general  form.  The  diagonal  elements  of  the  A  matrix  in 
the  Appendix  have  the  highest  probabilities  which  can  be 
interpreted  as  a  good  chance  of  recurrence  in  a  single  step. 
This  matrix  is  square  with  each  side  numbering  to  N,  the 
number  of  allowed  states.  The  number  of  observations  don't 
contribute  to  the  definition  of  this  matrix  because  it  deals 
strictly  with  the  states  of  the  model . 

The  initial  state  matrix  can  also  be  figured  by  beta 
density  function  integration.  In  the  case  of  the  3  state 
model,  starting  off  in  the  ambivalent  state  is  assumed  so 
the  beta  distribution  is  parameterized  to  allow  for  that. 
This  matrix  has  one  row  and  as  many  columns  as  allowed 


Markov  states.  The  row,  as  always,  has  to  sum  to  unity  to 
be  axiomatically  correct. 
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Observational  data 

The  aggregate  vector  fields  are  not  the  only  physical 
data  that  needs  to  be  collected.  Observing  the  bearing 
vectors  is  necessary  to  re-parameterize  the  HMM.  Here  it 
will  be  done  by  a  random  number  generator.  Appendix  B  shows 
the  MathCad  for  generating  observations.  The  number  of 
times  to  re-estimate  the  matrices  is  represented  by  the 
number  of  columns  in  the  observation  matrix  0.  The  re¬ 
estimation  algorithm  executes  after  each  set  of  data  is 
seen.  The  sets  of  data  correspond  to  the  columns  of  the  0 
matrix.  Thus,  by  controlling  how  many  columns  there  are  in 
the  matrix  and  the  length  of  the  columns,  we  can  decide  how 
often  to  re-estimate  the  model  parameters. 

By  asking  for  random  numbers  based  on  a  specific 
distribution  there  is  an  added  advantage.  This  is 
essentially  like  taking  the  observations  from  a  true  matrix 
biased  to  a  specified  state.  The  bias  state  can  be 
determined  by  how  we  parameterize  the  random  number 
generator.  Essentially,  we  will  see  if  the  formulation  is 
stable  later  on  by  giving  it  a  designed  matrix  and  seeing  if 
the  model  re-estimates  to  an  obvious  state.  If  the  most 
likely  state  is  similar  to  the  true  state  designated  by  the 
designed  data,  the  model  is  re-estimating  well.  These  tests 
are  included  in  the  results  chapter. 
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The  0  matrix  in  Appendix  B  can  be  interpreted  as  the 
data  that  the  satellite  gathers  as  information  before  a 
prediction  has  to  be  made  on  future  location.  It  can  be 
thought  of  as  data  that  is  from  the  distant  past .  For 
models  with  4  observations  and  different  numbers  of  states 
the  0  matrix  must  not  change.  Furthermore,  models  with  ,6 
observation  symbols  and  different  states  have  their  own  set 
of  data  because  the  definition  of  the  look  area  changes  and 
becomes  smaller.  This  was  done  because  one  data  set  could 
have  data  from  0-3,  while  the  other  data  set  could  range 
from  0-5.  The  sets  were  kept  the  same  across  the  state 
variations  so  each  formulation  could  be  re-estimated 
equally.  This  allows  for  good  comparisons  between  models 
later. 

The  data  that  will  be  used  to  do  the  comparing  is  the 
realO  matrix  also  in  Appendix  B.  This  data  can  be  thought 
of  as  being  in  the  recent  past.  This  matrix  also  varies 
between  models  with  different  observations  symbols.  Like 
the  0  matrix,  however,  it  is  the  same  across  the  number  of 
allowed  states.  Once  the  models  have  been  parameterized 
with  the  0  matrix  to  their  final  X  parameters,  the  realO 
matrix  will  be  used  in  comparing  the  models  with  different 
states. The  probability  of  the  entire  realO  matrix  will 
be  assessed  with  the  final  X  parameters.  The  model  that 

^he  final  matrices  can  be  thought  of  as  the  re-estimated  matrices  used  at  the  time  when  a  prediction  is 
needed. 

^^Remember  that  X  denotes  the  matrices  A,  B,  and  n  in  the  HMM  formulation. 
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predicts  this  realO  with  the  highest  probability  will  be 
considered  the  proper  model  to  minimize  modeling  error.  The 
number  of  columns  and  rows  can  change  from  the  realO  to  the 
0  matrix  to  demonstrate  some  versatility  in  the  model 
formulation.  The  comparison  between  order  of  models  can  be 
seen  as  part  of  the  results  chapter. 

Forward  and  Backward  Probabilities 

In  order  to  solve  for  the  probability  of  certain 
sequences  and  all  the  other  problems  associated  with  HMM, 
the  forward  and  backward  probabilities  have  to  be  solved. 
These  variables  come  from  the  fact  that  a  complete 
enumeration  of  probabilities  for  each  sequence  of  outcomes 
is  computational  infeasible.  The  calculation  is  on  the 
order  of  2T*n'^  [14]  ,  where  T  is  time  observed,  while  N  is 
the  number  of  states .  With  larger  values  of  N  and  T  the 
problem  quickly  expands  exponentially. 

The  forward  variable,  a^{±) ,  is  defined  as  the 
probability  of  a  partial  observation  sequence,  0^,  02,...,0t 
and  state  at  time  t,  given  the  model  X.  The  key  to  this 
concept  is  that  since  there  are  only  a  certain  number  of 
states,  regardless  of  how  long  the  sequence  is,  the  object 
is  constrained  to  transition  among  the  allowed  number  of 
states  throughout  the  sequence.  The  definition  the  forward 
variable  is : 

a,ii)  =  P{0„0„...,0,,q,=S,\A) 


(14) 
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To  solve  for  the  forward  variable,  two  steps  have  to  be 
taken.  The  first  step  includes  finding  the  joint 
probability  of  being  in  state  and  finding  observation  Oi. 

a,(0=^A(Oi)  (15) 

The  i  index  is  iterated  from  1  to  N.  The  second  step  is 
represented  by  Figure  3 .  It  shows  how  to  get  to  state  Sj 
from  Si  at  any  time  t+1.  This  step  inductively  finds  every 
forward  variable  after  the  initial . 


«,+iO)  = 


N 


L  »=i 


(16) 


The  j  index  is  iterated  from  1  to  N,  while  the  t  index  is 
iterated  from  l  to  T-1.  Summing  at(i)aij  over  every  state 
yields  the  probability  of  being  in  state  Sj  at  time  t+1. 

This  is  true  because  a^-aij  is  the  joint  probability  of 
seeing  Oi,  O2,...,  O-p  and  Sj  at  time  t+1  from  state  S^  at 
time  t.  The  calculation  effort  here  is  on  the  order  of  N^t. 
This  saves  orders  of  magnitude  in  the  number  of  calculations 
required  [14] . 

The  backward  variable  is  essentially  the  mirror  image 
of  the  forward  variable.  It  is  the  probability  of  the 
partial  sequence  of  observations  from  t+1  to  T,  given  state 
Si,  at  time  t  and  model  X.  The  backward  variable  can  be 
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Figure  3 .  Representation  of  Computations  required  for 

Forward  Variable  (Extracted  from  [14] ,  pg.  262) 
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given  by: 

A(0  =  P(0,„,0„,  ,...,0,\q,  =  5,  ,A) 


(17) 


Since  it  works  backward,  the  initialization  step  has  to  set 
the  equal  to  unity  for  states  1  to  N  as  in  equation 


(18)  . 


y9,(/)  =  l 


(18) 


The  induction  step  can  be  seen  in  Figure  4 .  It  shows  that 
in  order  to  be  at  state  Sj  at  time  t  and  to  account  for  the 
observation  sequence  from  time  t+1  to  T  we  have  to  account 
for  being  in  all  states  Sj  at  time  t.  The  induction  step 
can  be  seen  as  equation  (19) . 


A,(o=ZV)(o,*,)A*,a) 

^■=1  (19) 

The  aij  term  above  accounts  for  the  transition  from  Si  to 
Sj .  The  bj (Ot+i)  term  accounts  for  the  observation  seen  at 
time  t+1  in  state  j.  Finally,  the  Pt+l(3)  term  accounts  for 
the  rest  of  the  observation  sequence  from  state  j .  The  time 
step  starts  at  T-l  and  works  to  1  for  all  states  from  1  to 
N.  The  computational  effort  is  also  simplified  here  along 
the  same  lines  as  above  in  the  forward  variable  [14] . 

The  forward  and  backward  variable  are  not  numerically 
simple  to  compute.  It  is  not  difficult  to  see  that  for 
large  values  of  time  the  forward  and  backward  variables  will 
result  in  underflow  of  any  computer.  Even  the  small  scale 
model  shown  in  the  methodology  results  in  very  small 
probabilities.  Of  course,  the  underflowed  values  will  be  at 
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T  for  the  forward  variable  while  the  problem  for  the 
backward  variable  arises  at  the  beginning  of  time. 

The  actual  matrices  for  this  model  can  be  seen  as  part 
of  Appendix  C  under  FwdProb  and  BackProb,  respectively.  The 
solution  to  this  problem  comes  from  scaling  the  variables  to 
stay  within  the  dynamic  range  of  the  computer.  Essentially, 
the  scale  is  the  inverse  of  the  forward  variable  summed  for 
the  number  of  states  in  the  model  [14] ,  [16]  . 

The  scaling  factors  were  calculated  by  two  ways  in 
Appendix  C.  The  first  way  was  done  in  traditional  MathCad 
symbology  as  the  (scale)  variable.  The  second  way  was 
accomplished  using  the  programming  structures  within  MathCad 
and  can  be  seen  as  simply  (c)  in  Appendix  C.  Both  ways  were 
done  to  assure  the  procedure  was  correct.  Improper  scaling 
could  seriously  flaw  the  reliability  of  final  results. 

Once  we  were  assured  that  the  numbers  were  correct, 
they  could  be  used  in  calculations.  By  multiplying  the 
scale  by  each  forward  and  backward  variable  throughout 
calculations  they  will  cancel  each  other  out  by  divisions. 

It  is  a  convenient  way  to  deal  with  the  underflow  problem 
and  can  help  us  to  calculate  sequence  probabilities  later  on 
[14]  ,  [16]  . 


Most  likely  State 

Knowing  what  Markov  state  the  model  is  in  at  a 
particular  time  is  difficult  because  it  is  not  seen.  To 
recall  the  bias  coin  example,  it  would  be  like  finding  out 
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what  bias  coin  is  being  used  based  on  only  a  set  of  heads 
and  tails.  For  the  tracking  problem  being  discussed  here, 
the  bearing  observations  will  provide  the  only  clue  as  to 
what  protocols  the  object  is  using. 

The  mathematics  behind  determining  this  can  be  very 
complex.  The  Viterbi  algorithm  is  one  way  to  find  the  best 
sequence  of  states  for  any  given  set  of  observations.  It 
maximizes  the  probability  of  the  single  state  sequence  given 
the  observations  and  the  parameters  of  the  model  (max  P(Q|0, 
X)  [14] ,  [17]  .  Although  there  are  plenty  of  HMM  applications 
that  use  the  above  algorithm,  the  most  probable  sequence  is 
of  little  use  in  this  application.  The  current  state  is 
more  important  for  this  situation.  That  way  trends  for 
future  movement  can  be  verified  against  what  the  believed 
protocols  are  at  that  time. 

In  order  to  find  out  the  most  likely  state  at  each 
time,  we  want  to  maximize  the  expected  number  of  correct 
individual  states  over  a  time  period.  To  do  this,  we  need 
to  find  the  probability  of  being  in  state  Sj[  at  time  t, 
given  the  sequence  0  and  the  model  X. 

The  formal  definition  is  [14] : 

r,ii)  =  Pi^,=S,\0,A)  (20) 

Appendix  C  shows  how  it  can  be  defined  in  terms  of  forward 
and  backward  variables.  It  is  a  probability  measure  so 
summing  over  states  for  each  time  yields  unity.  Verifying 
this  in  the  calculations  indicates  the  mathematics  are  sound 
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up  to  this  point.  Errors  in  definition  and  indices  would 
yield  results  other  than  unity  for  the  summation. 

The  measure  serves  to  show  the  most  likely  state  at  a 
particular  time.  The  problem  with  the  calculation  is  that 
it  doesn't  have  any  regard  for  the  probability  of  occurrence 
of  sequences  of  states.  In  other  words,  it  could  solve  for 
a  bogus  sequence  of  states  that  aren't  even  possible  because 
the  state  transitions  aren't  feasibly  allowed  [14].  As 
mentioned  above,  the  current  state  is  the  most  important 
information,  while  the  sequence  is  of  little  use.  This 
allows  us  to  use  the  criteria  for  most  likely  expected  state 
at  individual  times  rather  than  expected  state  sequence  of 
the  Viterbi  algorithm. 

As  an  additional  check  on  the  scaling ‘factors,  see  both 
calculations  of  y  in  Appendix  C.  With  and  without  the 
scaling  factors  the  matrices  should  be  equal.  Since  they 
both  are,  the  calculations  are  correct.  To  actually 
determine  the  most  likely  state,  look  down  each  column 
(time)  and  across  each  row  (Markov  state)  in  Appendix  C. 

For  each  time  (column)  there  will  be  a  row  (state)  with  the 
highest  probability.  This  will  be  the  state  that  the  object 
is  expected  in  and  the  number  will  show  what  that 
corresponding  probability  is.  Obviously,  if  there  is  not  a 
state  that  has  a  high  probability  it  will  be  hard  to  draw 
useful  information  from  the  y  calculation  in  equation  (20)  . 


38 


Re-estimation  of  Parameters 

The  re-estimation  of  parameters  is  the  most  tedious  and 
difficult  problem  to  solve  in  the  study.  It  is,  however, 
the  most  important.  Without  being  able  to  modify  the 
parameters  of  the  model,  observations  would  be  useless.  Any 
predictions  of  movement  would  have  to  be  based  on  the 
initial  parameters  estimated  before  any  data  was  taken.** 

This  defeats  the  purpose  of  observing  in  the  first  place. 

The  algorithm  for  re-estimating  can  be  seen  as  part  of 
Appendix  C  in  the  matrix  routine.  Yet,  another  variable  can 
be  defined  to  help  validate  the  formulation.  The 
probability  of  being  in  state  Si  at  time  t,  and  state  Sj  at 
time  t+1,  given  the  model  and  observation  sequence  is 
defined  as: 

Writing  this  equation  in  terms  of  forward  and  backward 
variables  is  part  of  Appendix  C.  Again,  if  all  formulation 
and  mathematics  work  has  been  done  correctly  up  to  this 
point  then: 


(22) 

For  the  time  0  of  the  model  shown  in  the  calculations  of 
Appendix  C,  the  left  hand  and  right  hand  side  of  equation 
(22)  equal  each  other.  A  slight  change  of  notation  occurs 


**Remember  that  the  parameters  of  the  model  (i)  are  given  Ity  the  matrices  A,  B,  tc. 
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here.  The  ^  of  equation  (21)  is  represented  by  a  Sat  zero  in 
the  appendix.  It  will  become  clear  as  to  why  later  on. 

The  matrix  algorithm  looks  more  complicated  than  it  is. 
All  the  small  routines  are  set  up  within  one  large  while 
loop.  This  while  loop  is  conditioned  on  being  less  than  K 
to  operate.  Recall  from  Appendix  B  that  K  is  defined  as  the 
number  of  columns  in  the  observational  data  (0  matrix) . 

Also,  recall  that  the  observations  are  set  up  in  a  matrix  so 
a  re-estimation  can  take  place  after  each  column.  The  while 
loop  here  is  what  causes  a  re-estimation  every  column. 

Inside  the  main  loop  several  variables  are  initialized. 
The  next  major  step  in  the  algorithm  is  the  calculation  of 
S.  The  definition  of  S  is; 

^=1  (23) 

This  is  the  expected  number  of  transitions  from  to  Sj 
[14] .  S  has  time  numbered  columns  and  state  numbered  rows. 
This  calculation  is  a  triple  nested  loop  for  every  element 
of  the  matrix  and  every  time.^^  Trying  to  build  the  ^  matrix 
of  equation  (21)  is  difficult  because  it  has  essentially 
three  dimensions.  This  poses  an  implementation  problem  for 
MathCad  because  it  will  not  recognize  such  a  matrix  as 
valid.  Therefore,  the  E  matrix  had  to  be  constructed  at 
each  iteration  of  the  E  loop.  The  sumZeta  notation  simply 
handles  adding  the  appropriate  time  elements  together. 

^^The  above  Satzero  was  defined  that  way  because  it  was  at  the  first  time  possible.  When  dealing  with 
only  one  time  the  summation  can  disappear  and  the  S  and  ^  variables  equal  each  other. 
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Notice  how  the  variables  for  each  separate  loop  are  re¬ 
initialized  at  the  end  of  the  loop  rather  than  at  the 
beginning . 

Once  the  S  matrix  has  been  constructed,  the  y  matrix 
can  be  calculated  for  the  while  loop.  This  is  the  same 
variable  as  defined  before  in  equation  (20)  .  Now  that  we 
are  in  a  loop,  this  value  has  to  be  kept  current  with  the  X 
at  the  corresponding  time.  The  previous  work  was  for 
purposes  of  seeing  what  was  going  on  qualitatively  and 
making  sure  the  mathematics  were  constructed  properly.  The 
y  matrix  will  be  indirectly  calculated  from  X.  Therefore, 
using  the  X  from  previous  work  (i.e.,  outside  the  loop) 
would  yield  incorrect  results. 

With  the  two  above  matrices  constructed,  a  re-estimated 
Aij  matrix  can  be  determined.  The  operation  is  given  by: 

Sf.ov) 

A  =  Izi _ 

'=1  (24) 

The  new  A  matrix  requires  summing  over  i  and  j  as  seen  in 
equation  (24)  .  Essentially,  the  new  A  matrix  is  the 
expected  number  of  transitions  from  state  to  state  Sj 
divided  by  the  expected  number  of  transitions  from  state  Si 
[14]  . 
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The  new  7C  matrix  is  almost  trivial .  It  is  taken  simply 

as : 

^<  =  r,0)  (25) 

It  is  the  expected  frequency  in  state  Sj  at  time  t=l  [14] . 

The  new  B  matrix  is  much  more  complicated.  It  will 
require  the  operation: 

tr,0) 

"o  ^  sJ.O^=Vk 
T 

Zr,(') 

'=1  (26) 

This  variable  is  essentially  the  expected  number  of  times  in 
state  j  and  observing  symbol  V]^  divided  by  the  expected 
number  of  times  in  state  j  [14] .  First,  notice  the  index  on 
the  time.  It  is  from  T  rather  than  t-1.  This  is  because 
the  B  matrix  handles  observations  which  occur  at  every  time. 
The  last  state  transition  occurs  at  t-1.  Because  the  A 
matrix  handles  state  transitions,  its  index  can  only  go  as 
high  as  t-1.  Secondly,  notice  the  note  under  the  top 
summation.  This  summation  is  conditioned  on  the  symbol 
observed.  In  MathCad,  the  Y  variable  handles  this 
summation.  The  column  of  Y  to  get  the  y  added  to  it  is  the 
same  column  corresponding  to  the  observation  symbol  seen. 

The  observed  symbol,  of  course,  comes  from  the  training  data 
0  matrix.  The  element  depends  on  where  K  is  and  at  what 
time  the  loop  is  looking  at . 
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Another  implementation  issue  arises  in  this  algorithm, 
however.  Because  no  data  set  can  ever  be  complete,  the  re- 
estimation  algorithm  can  conceivably  create  elements  of  a  B 
matrix  that  are  equal  to  zero.  If  the  training  set  never 
sees  one  symbol  it  will  re-estimate  that  symbol  probability 
to  0 .  That  is  never  really  true  so  a  default  value  (s)  has 
to  be  given  to  that  matrix  element  in  the  absence  of 
training  data.  Once  that  is  done,  the  row  has  to  be 
rescaled  so  that  the  sum  of  the  row  is  equal  to  unity  [14] . 

Once  all  the  matrices  have  been  re-estimated,  they  have 
to  be  reset  to  the  old  matrices  so  the  loop  can  proceed  with 
the  new  parameters.  Once  the  final  column  of  training  data 
is  processed,  the  final  re-estimates  of  the  parameters 
become  the  output  matrices .  MathCad  will  only  output  one 
matrix  at  a  time  so  all  three  matrices  are  combined  for 
output.  Once  out,  the  aggregate  matrix  is  separated  into 
its  components  again. 

Of  course,  all  of  the  axioms  of  probability  must  be 
obeyed  to  have  meaningful  parameters  of  the  model.  All  of 
the  rows  in  each  of  the  matrices  have  to  be  equal  to  unity 
when  summed.  This  check  will  serve  as  a  test  of  the 
algorithm's  correctness.  It  has  been  proven  that  using  the 
procedure  in  Appendix  C  will  either  generate  new  A,  B,  and  n 
matrices  that  increase  maximum  likelihood  of  a  sequence  or 
stay  at  an  already  critical  point  in  the  likelihood  function 
[14]  . 
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Determining  the  best  observation  sequences 

Once  the  parameters  of  the  model  have  been  re¬ 
parameterized  finally,  the  highest  probability  observation 
sequences  may  be  determined.  The  probability  of  a  given 
observation  sequence  is  given  by: 

I  N 


P{0 


^)  =  Z«r(0 


(27) 


These  probabilities  can  get  extremely  small  for  long  time 
durations.  In  that  case,  using  log  probabilities  will  be 
best  because  computer  underflow  becomes  an  issue. 
Fortunately,  the  time  durations  used  don't  cause  these 
problems  so  the  calculation  can  be  direct. 

Time  durations  are  kept  small  because  of  the  necessity 
for  this  problem  to  come  up  with  an  answer  in  a  short  time. 
If  the  algorithm  took  days  and  days  of  observation  to 
analyze  the  object,  it  would  probably  not  be  any  better  than 
conventionally  searching  and  destroying  the  target. 

The  essence  of  this  new  problem  revolves  around  finding 
the  probabilities  of  the  most  likely  observation  sequences. 
Looking  at  equation  (27)  again,  the  forward  variable 
calculation  has  to  be  accomplished  first.  Appendix  D  starts 
out  with  this  exact  calculation.  Notice  that  the  scale  is 
not  included . 

After  that,  the  probability  of  each  sequence  has  to  be 
determined.  For  very  small  models  this  can  probably  be  done 
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manually  although  it  would  be  very  tedious.  For  any  model 
the  possible  number  of  combinations  of  symbols  is  given  by: 


(n)  denotes  the  number  of  possible  observations.  (r) 
denotes  the  length  of  the  observation  set  [18] .  For  this 
model  we  have  4  possible  observations  and  are  allowing  the 
set  to  be  3  symbols  long.  By  doing  the  calculation,  there 
are  64  different  sets  of  observations  to  extract  the 
probability  from.  For  models  of  more  complexity  this  number 
goes  up  exponentially. 

The  sequence  routine  in  Appendix  D  generates  all  the 
possible  sequences  for  our  model.  It  requires  as  many 
nested  loops  as  there  are  number  of  observations  allowed  in 
the  length  of  the  set .  Because  it  allows  the  same  number  of 
sequences  as  equation  (28)  the  procedure  was  done  correctly. 

The  LOGP  routine  does  the  actual  calculation  of 
probabilities.  It  uses  the  forward  probabilities  generated 
at  the  beginning  of  the  Appendix  to  output  a  probability  for 
each  sequence.  The  Bestset  routine  is  a  combination  of  both 
sequence  and  LOGP  routines  to  generate  a  probability  of  each 
observation  sequence .  Part  of  the  output  can  be  seen  below 
the  routine.  As  we  can  see,  the  code  starts  to  enumerate 
all  the  sequences  and  find  the  probability  associated  with 
each  above  it.  The  0,0  element  shows  how  many  total 
sequences  there  are  and  the  0,1  element  shows  what  the  final 
summed  probability  across  row  l  is. 
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This  summed  probability  should  sum  to  one.  This  is 
because  we  are  looking  at  the  entire  event  space.  The 
reason  why  the  sum  isn't  unity  is  numerical.  Allowing  the 
time  to  only  go  to  2  iterations  will  indeed  sum  the  row  to 
unity.  Allowing  the  time  to  go  to  4  iterations,  however, 
will  further  degrade  the  summation.  After  much  deliberation 
it  was  determined  that  because  of  the  small  nature  of  the 
probabilities  and  the  sheer  number  of  them  needing  to  be 
added  together,  it  is  entirely  possible  to  generate  the 
error  MathCad  shows.  This  conclusion  is  consistent  with 
allowing  the  time  to  vary  and  observing  how  the  summation 
gets  closer  to  and  farther  from  unity. 

For  all  practical  purposes,  the  end  result  will  not  be 
affected.  Relatively  speaking,  the  sequence  order  will  not 
change  because  all  probabilities  will  be  affected  by  the 
same  issue.  One  way  to  deal  with  it,  therefore,  is  to 
simply  rescale  all  the  probabilities.  Dividing  all 
probabilities  by  the  sum  of  probabilities  will  yield  numbers 
that  will,  in  fact,  sum  to  unity.  The  Bestset  routine 
contains  such  a  normalizing  operation  so  that  the  final 
probabilities  do  sum  to  unity. 

Once  all  the  sequences  have  been  identified  and  given  a 
probability  with  Bestset,  they  are  sorted  out  to  find  the 
most  likely  sequences.  The  rsort  function  in  MathCad  does 
that  for  us.  It  sorts  the  matrix  from  smallest  to  largest 
based  on  the  first  row  of  the  column.  Taking  the  transpose, 
reversing  it,  and  taking  another  transpose  was  the  easiest 
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way  to  get  the  program  to  sort  from  largest  to  smallest. 
MathCad  doesn't  have  a  function  that  will  sort  from  largest 
to  smallest  automatically.  The  rsort  routine  is  based  on  a 
heapsort  algorithm  found  in  Numerical  Recipes  in  C  [19] . 

The  new  matrix  can  be  seen  below  the  manipulations .  It 
sorts  the  most  likely  sequences  from  most  probable  to  least 
probable . 

Determining  the  best  model  and  seeing  trends 

Once  the  most  likely  sequences  have  been  sorted  out, 
trends  can  be  determined  in  how  the  object  will  move.  An 
incorrect  model  will  lead  to  poor  information  about  the 
object,  however.  On  the  other  hand,  knowing  the  proper 
order  will  improve  the  quality  of  the  information.  This 
section  is  based  primarily  on  trying  to  find  the  best  model 
given  the  limited  training  sequences  used  to  re-parameterize 
the  model . 

Up  to  this  point,  all  observations  have  been  from  the 
observed  data  in  the  0  matrix.  This  is  distant  past  data 
that  is  considered  given.  The  realO  matrix  is  data  that  was 
seen  in  the  recent  past  (see  Appendix  B) .  This  division  of 
the  data  is  arbitrary  and  used  to  help  verify  the  model 
order.  It  is  the  realO  observations  that  will  be  used  to 
determine  what  model  is  most  appropriate  for  that  data.  It 
is  important  to  realize  that  we  should  not  use  the  same  0 
matrix  to  determine  the  order.  The  0  matrix  was  used  to 
parameterize  the  model  so  the  parameters  reflect  that  matrix 
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specifically.  A  new  set  of  data  has  to  be  given  to  show 
best  order  most  appropriately. 

Unfortunately,  when  the  number  of  observation  symbols 
for  a  model  changes,  this  matrix  will  have  to  be 
regenerated.  This  is  because  one  model  won't  be  able  to 
handle  the  same  number  of  symbols  that  another  model  might 
demand.  The  number  of  observation  symbols  for  this 
application  of  HMM  can  be  tailored  to  fit  how  much  data  is 
available  and  to  what  degree  of  accuracy  we  need.  Computer 
speed  limited  this  study  to  only  4  and  6  symbols.  If  it 
were  deemed  necessary,  however,  every  degree  on  the  compass 
could  have  its  own  symbol. 

Throughout  the  state  variations,  however,  the  realO 
matrix  will  remain  unchanged.  The  state  is  an  important 
area  to  look  at  when  determining  order  because  this  is  the 
part  of  the  model  that  can't  be  set  arbitrarily. 

Furthermore,  the  state  is  also  the  part  of  a  HMM  that  can't 
be  observed.  That  makes  it  impossible  to  determine  an  order 
based  strictly  on  seeing  what  the  output  is.  Therefore,  we 
have  to  find  a  way  to  determine  an  order  based  on  the 
limited  observational  data  from  the  realO  matrix. 

The  realO  matrix  and  calculations  for  the  order  can  be 
seen  as  part  of  Appendix  E.  The  premise  is  as  follows:  The 
probability  of  each  sequence  of  realO  is  determined  by  LOGP 
routine  in  previous  section.  The  LOGP  routine  is  based  on 
the  parameters  of  the  model  in  question.  If  these 
parameters  are  closer  to  ideal  than  other  models,  higher 
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probabilities  will  be  determined  for  the  sequences  in  realO. 
This  is  true  because  more  ideal  parameters  would  predict  the 
sequences  in  realO  with  a  higher  probability.  If  we  were  to 
take  the  log  of  these  parameters,  the  above  strategy  changes 
slightly.  The  highest  probability  will  now  be  closest  to 
zero.  A  set  of  the  highest  probabilities,  when  summed,  will 
be  closer  to  zero  than  another  set  of  probabilities  summed. 

In  essence,  a  score  is  generated  to  evaluate  a  model. 
The  score  closest  to  zero  on  the  log  scale  will  determine 
the  best  model.  Remember,  this  is  true  because  the 
probabilities  come  straight  from  the  parameters  of  the 
model.  It  is  not  unlike  doing  a  reduced  chi  squared  fit. 

We  try  to  minimize  the  vertical  distance  between  the  data 
and  line  drawn.  The  summed  vertical  distance  is  the  reduced 
chi  square  and  once  that  is  minimized,  the  best  fit  is  known 
[20] .  Our  line  of  fit  is  set  at  the  x-axis.  The  Markov 
states  of  the  model  are  the  only  things  altered  to  try  to 
fit  that  line. 

It  is  important  to  note  that  this  isn't  really  a  chi- 
squared  fit.  It  is  only  like  a  chi-squared  fit.  Chi- 
squared  has  slopes  and  uncertainties  that  play  a  role, 
whereas  the  vertical  separation  is  the  only  important  part 
of  the  above  procedure.  Also,  this  model  doesn't  move  the 
line  to  reduce  the  score,  the  model  is  fundamentally  changed 
so  the  data  will  fit  the  line  better.  It  is  the  inverse 


process . 
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Changing  the  states  of  the  model,  however,  requires 
that  the  training  data  not  change.  If  the  data  were  to 
change  then  comparing  two  different  models  wouldn't  have  any 
significance.  That  is  the  reason  realO  was  kept  constant 
for  a  given  set  of  allowed  observation  symbols.  Finally, 
after  all  the  new  calculations  the  different  states  can  be 
compared  based  on  the  realO  data.  The  lowest  score  (i.e., 
lowest  summed  vertical  difference)  will  best  determine  the 


order  of  the  model . 


CHAPTER  IV 


MODEL  FORMULATIONS 

As  mentioned  above,  there  has  only  been  one  formulation 
of  the  model  discussed.  We  have  been  using  3  Markov  states 
and  4  observation  symbols  to  explain  the  algorithm  in  the 
previous  chapter.  The  3  states  correspond  to  a  speed  state, 
conceal  state,  and  ambivalent  state.  Basically,  the  Markov 
states  represent  the  mindframe  that  the  tracked  object  is 
in.  The  4  observation  symbols  correspond  to  each  direction 
of  the  compass.  This  made  the  methodology  simpler  to 
explain  because  there  was  a  physical  significance  associated 
with  each  part  of  the  model.  For  this  application,  perhaps 
there  are  more  sophisticated  states  that  we  aren't 
considering.  Also,  there  might  be  more  observation  symbols 
than  we  initially  allowed. 

Under  these  circumstances,  the  model  discussed  up  to 
now  would  become  useless.  Modeling  error  would  be  large 
[10] .  It  wouldn't  be  an  accurate  estimate  of  what  is 
happening.  The  purpose  of  this  chapter  is  to  try  to  go  past 
what  we  can  lend  physical  significance  to  in  Markov  state 
modeling  and  determine  the  best  order  of  the  model . 
Observation  symbol  variation  will  also  be  explained. 
Regardless  of  the  change,  however,  the  matrices  should  stay 
axiomatically  correct.  All  rows  still  have  to  sum  to  unity 
in  order  to  be  valid  when  used  later  on. 
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Only  the  alterations  will  be  discussed  below.  It  would 
be  redundant  to  discuss  every  model  totally  since  there  is 
much  overlap  between  the  models.  The  interesting  parts  of 
each  model  are  how  they  differ  from  the  fundamental  model 
used  in  the  methodology  chapter.  A  list  of  each  model  can 
be  seen  as  part  of  table  1 . 


MODEL  NUMBER 

NUMBER  OF 

STATES 

NUMBER  OF  OBSERVATION 

SYMBOL 

1 

3 

4 

2 

3 

6 

3 

4 

4 

4 

4 

6 

5 

5 

4 

6 

5 

6 

Table  1:  List  of  each  Model  Formulation  used. 


Six  Observations 

In  this  study,  4  and  6  observations  were  examined.  For 
this  application,  the  number  of  symbols  is  really  only 
limited  by  the  computer  speed.  The  issue  is  the  precision 
needed  to  determine  where  the  object  will  be  at  a  specific 
time.  Theoretically,  each  degree  of  the  compass  can  be  used 
if  needed,  but  realistically  that  doesn't  make  any  sense. 

A  trend  in  movement  will  be  predicted  in  as  few 
observations  as  possible.  This  study  will  use  3 
observations.  Many  observation  symbols  requires  more  data. 
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With  a  very  complicated  model,  the  time  it  takes  to  predict 
a  movement  might  be  longer  than  just  going  out  and  looking 
for  it  with  no  prior  knowledge  at  all.  With  that  in  mind, 
it  might  not  be  very  useful  to  go  much  more  complicated  than 
demonstrated. 

The  changes  needed  to  increase  the  number  of 
observation  symbols  can  be  seen  as  part  of  Appendix  F.  The 
convenient  aspect  of  this  model  is  that  it  is  set  up  in  way 
that  takes  whatever  matrices  given  to  it  and  does  the 
complicated  work  without  being  altered.  In  other  words,  the 
initial  matrices  section  of  the  methodology  is  the  only  part 
of  the  model  that  needs  to  be  changed  to  change  the  order. 
The  rest  of  the  methodology  is  general  enough  to  adopt 
itself  to  given  matrices. 

The  same  vector  field  is  used  in  the  6  symbol  and  4 
symbol  models.  The  major  change  happens  in  the  number  of 
directions  looked  in.  The  encompassed  angle  is  decreased 
from  90  to  60  degrees.  In  accordance  with  that,  fewer  rows 
or  columns  of  the  field  are  examined  for  any  one  direction. 
Also,  the  weighting  had  to  be  changed  slightly  to  account 
for  this  decrease. 

Other  than  that,  the  probabilities  for  each  of  the 
states  are  calculated  the  same  way  as  in  the  4  symbol  model. 
The  observation  symbol  matrix,  B,  increases  columns  from  4 
to  6  to  fit  the  new  symbol  probabilities.  The  A  and  tc 
matrices,  however,  remain  untouched  because  the  number  of 
Markov  states  doesn't  change  at  all.  Changing  the  number  of 
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symbols  is  really  not  that  difficult.  The  rest  of  the 
formulation  will  pick  up  these  matrices  and  operate 
according  to  how  the  methodology  stated  in  the  previous 
chapter . 


Four  States 

Changing  the  number  of  allowed  Markov  states  is 
unfortunately  not  as  trivial.  As  mentioned  above,  we  have 
only  talked  about  a  3  state  model .  Each  one  of  these  had  a 
physical  significance.  Here  we  have  4  states.  There  are  4 
symbols  so  as  to  not  confuse  the  issue  more  than  it  needs  to 
be.  Now  that  we  know  the  formulation  is  correct  (see  result 
chapter)  physical  significance  can  play  less  of  a  role  in 
the  Markov  states .  The  changes  made  can  be  seen  as  part  of 
Appendix  G. 

From  the  given  vector  field,  the  4  directions  are 
calculated  and  applied  to  the  speed  and  conceal  state.  The 
speed  and  conceal  states  were  kept  in  this  formulation  as  a 
creative  choice.  The  ambivalent  state,  however,  goes  away. 
Instead,  we  parameterize  a  beta  probability  density  function 
and  split  the  ambivalent  state  into  two  separate  states. 

The  parameters  are  also  a  creative  choice.  The  parameters 
were  designed  to  make  the  two  states  compliment  each  other. 
What  is  the  physical  meaning  to  these  2  states?  There  might 
not  be  any.  It  is  not  important.  Meaning  can  be  examined 
and  determined  later  if  this  formulation  is  the  best  model . 
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Parameterizing  these  states  with  observational  data  will 
adjust  symbol  probabilities  accordingly. 

What  about  state  transition  probabilities  and  initial 
state  probabilities?  Even  though  there  are  more  states  than 
before,  this  estimation  will  not  change.  State  transitions 
will  be  approximated  using  beta  density  functions  set  up  to 
be  recurrent  in  one  iteration.  The  initial  state 
probabilities  will  also  be  estimated  with  a  beta  function 
parameterized  to  be  more  probable  for  the  semi -ambivalent 
states . 

Unlike  the  symbol  changes  mentioned  above,  changing  the 
states  changes  all  the  matrices.  The  B  matrix  will  take  on 
another  row  to  accommodate  the  extra  state .  The  A  matrix 
will  be  4x4  now  instead  of  3x3 .  This  is  so  that  it  can  show 
the  extra  state  transitions .  The  n  matrix  will  take  on 
another  element  to  show  an  initial  probability  for  the  extra 
Markov  state . 


Five  States 

Here,  instead  of  4  Markov  states,  the  model  is 
formulated  with  five  states.  The  only  difference  comes  from 
how  we  formulate  the  states.  The  beta  function  is  used 
again  with  the  same  parameters  but  a  new  state  has  entered 
the  picture.  Here,  we  chose  to  bring  the  true  ambivalent 
state  back.  The  extra  state  uses  a  uniform  distribution  as 
the  first  estimate.  The  alteration  for  5  states  can  be  seen 
in  Appendix  H. 


CHAPTER  V 


RESULTS 

As  mentioned  above,  the  principle  behind  all  this  work 
has  been  to  come  up  with  a  novel  approach  to  determine  the 
order  of  a  HMM  so  as  to  predict  future  movements.  With 
those  predictions  we  can  then  hopefully  reduce  a  search  area 
when  trying  to  destroy  an  object.  The  model  has  to  be 
validated,  however,  before  any  useful  information  can  be 
taken  from  it 


Correctness  of  Model 

Essentially,  this  model  is  for  forecasting.  How  do  we 
go  about  validating  a  forecasting  model?  A  perfectly  good 
approach  is  to  take  data  from  the  distant  past  and  try  to 
predict  the  data  that  is  in  the  recent  past  [9] . 

The  simplest  way  to  test  the  model  is  to  design  several 
simple  sequences  to  force  a  certain  state  prediction.  If 
the  model  predicts  the  correct  state  based  on  the  sequence, 
it  is  stable.  We  can  design  trivial  secfuences  to  force  a 
prediction  to  each  state  of  the  fundamental  model.  The 
trivial  sequences  represent  the  distant  past,  while  the 

t 

state  being  tested  represents  the  recent  past.  The 
fundamental  model  was  used  to  test  the  formulation.  We 
don't  need  to  check  the  rest  of  the  models  because  they  are 
just  derivations  on  the  first.  Remember  that  all  other 
models  only  change  how  the  initial  matrices  are  estimated. 


56 


As  long  as  these  matrices  sum  rows  to  unity  and  the 
fundamental  model  works,  all  of  the  rest  are  correct  also. 

Table  2  has  the  sequences  chosen. It  should  be 
obvious  by  looking  at  Figure  2  or  Appendix  A  that  the  first 
sequence  should  force  predictions  to  the  conceal  state. 
Remember  that  the  highest  parts  of  the  field  correspond  to 
areas  of  more  concealment.  The  second  sequence,  on  the 
other  hand,  stays  with  the  lowest  parts  and  should  force  a 
speed  state.  The  lowest  parts  of  the  field  corresponds  to 
fast  areas  such  as  roads  and  clear  weather.  The  third 
sequence  makes  a  circle  and  obviously  forces  an  ambivalent 
state . 


Forced 

State 

Conceal 

Speed 

Ambivalent 

Observ. 

1 

0 

3 

0 

Observ. 

2 

0 

3 

1 

Observ. 

3 

0 

0 

2 

Table  2 :  Sequences  to  Determine  Correctness  of  Model 


Appendix  I  has  the  results  for  this  test.  By  examining 
the  A,  B,  and  n  matrices  we  can  see  where  the  model  predicts 
future  movements  and  state  transitions  based  on  highest 
probabilities,  y,  for  example,  shows  what  state  is  expected 
at  each  time.  The  highest  probability  should  be  with  forced 
state.  Also,  we  can  look  at  what  is  predicted  for  the 


^^Remember  that  0  is  north,  1  is  east,  2  is  south,  and  3  is  west 


57 


future.  If  the  most  probable  movement  predictions  mimic  our 
trivial  sequence  then  we  know  the  model  works .  Based  on  the 
results,  it  is  safe  to  say  that  this  model  does  force  the 
appropriate  states.  Hence,  the  formulation  is  stable  and 
correct . 

The  format  of  Appendix  I  is  also  the  format  that  we 
would  use  when  trying  to  attain  information  about  any 
observations.  The  A,  B,  and  7i  matrices,  y,  and  future 
predictions  matrix  should  provide  an  analyst  with  enough 
information  about  a  target  to  make  decisions  on  its  future 
movements.  In  addition,  the  three  of  these  results  when 
taken  together  should  also  be  enough  for  someone  to  make  a 
decision  on  whether  or  not  the  data  is  inconclusive. 

Order  of  Model 

Now  that  we  are  confident  that  the  basic  algorithm  is 
correct,  we  can  look  at  determining  the  proper  order.  Much 
deliberation  went  into  finding  a  method  that  would  show  the 
best  order  of  a  model  given  the  limited  training  data.  As 
mentioned  before,  an  obvious  way  is  to  define  a  truth  model 
and  generate  data  from  that.  Then,  see  if  the  model  could 
parameterize  to  that  truth  model  over  time.  There  are 
several  ways  to  test  the  model  depending  on  the  application 
[21,5]  . 

The  problem  arises  in  the  limited  amount  of  data.  A 
truth  model  could  forseeably  give  misleading  data  over  short 
time  spans.  This  would,  in  turn,  cause  the  model  to 
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parameterize  poorly  and  make  us  think  that  the  formulation 
is  wrong  when  it  isn't.  Therefore,  we  need  something 
different. 

Instead  of  a  truth  model,  consider  an  ideal  model. 

This  ideal  model  will  never  mislead  over  any  time  span.  The 
parameters  of  the  ideal  model  are  set  up  to  perfectly 
predict  the  data  that  was  given.  In  other  words,  plotting 
the  log  probabilities  of  the  sequences  would  be  very  close 
to  the  y=0  line  for  the  ideal  model .  The  summed  vertical 
distances  of  each  point  from  the  x-axis  would  be  very  close 
to  zero. 

Our  challenge  simply  becomes  finding  the  same  log 
probabilities  for  the  generated  models.  The  model  with  the 
smallest  summed  vertical  distances  is  closest  to  the  ideal 
model.  In  other  words,  that  model  will  have  the  most 
appropriate  order. 

Appendix  J  has  the  results  of  both  the  4  observation 
models  and  the  6  observation  models.  The  matrix  for  each 
model  contains  the  sorted  log  probabilities  of  the  sequences 
for  the  realO  matrix.  These  log  probabilities  are  based  on 
the  final  parameters  of  the  corresponding  model.  Remember 
that  the  realO  matrix  will  change  between  4  and  6 
observations  to  allow  for  the  possibility  of  more  symbols. 

Each  matrix  is  then  summed  to  produce  a  "score"  for 
that  model.  The  model  that  is  closest  to  the  ideal  will 
have  a  score  closest  to  zero  because  the  correct  model  would 
have  predicted  the  realO  matrix  perfectly.  The  score  for  a 
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set  of  observations  can  then  be  graphed  versus  the  state  of 
the  model.  Figure  5  shows  the  4  symbol  score,  while  Figure 
6  shows  the  6  symbol  scores . 


Model  Score  versus  Allowed  State 


Figure  5 .  4  Observational  Symbol  Scores  versus  Markov  State 


Model  Score  versus  Allowed  State 


Figure  6 .  6  Observational  Symbol  Scores  versus  Markov  State 
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As  the  number  of  states  increased,  the  score  was  better 
for  the  4  symbol  model.  This  would  indicate  that  5  states 
best  describes  the  realO  data.  As  the  states  increased  for 
the  6  symbol  model,  however,  there  was  a  plateau  effect. 

For  all  intensive  purposes  one  could  make  an  argument  that 
the  4  and  5  state  models  acted  the  same.  Both  were  better 
than  the  3  state  model,  however.  More  testing  would  have  to 
be  done  to  reach  a  conclusion  on  the  order  of  the  model  in 
both  cases. 

The  4  symbol  model  had  lower  scores  than  the  6  symbol 
model.  This  is  why  we  couldn't  compare  the  two  directly. 

As  we  mentioned  before,  more  observations  possible  means 
inherently  lower  probabilities.  Without  thinking  carefully 
about  it,  we  might  have  compared  the  two  and  decided  that 
the  4  symbol  models  were  automatically  better.  This  isn't 
necessarily  the  case.  The  nature  of  the  two  are  totally 
different  in  that  the  more  symbols  we  allow  the  more  data  we 
will  need.  This  is  a  definite  trade-off  to  consider  when 
deciding  how  many  symbols  is  appropriate  for  a  specific 
application. 

It  is  unfortunate  that  further  testing  couldn't  be 
conducted  to  reach  more  of  a  conclusion.  The  barrier  was  in 
the  run  time  it  took  for  each  model .  The  differences 
between  simple  and  complex  were  marked.  The  fundamental  3 
state,  4  observation  model  took  approximately  18  minutes  to 
run.  The  most  complex  model  with  5  states  and  6 
observations  took  on  the  order  of  8  hours  to  run.  The 


number  of  states  was  the  main  factor  in  anin  time,  while 
number  of  observation  symbols  was  inconsequential  in 
determining  how  long  the  model  would  run. 


CHAPTER  VI 


FURTHER  RESEARCH 

This  paper  is  by  no  means  all  inclusive.  There  are 
many  areas  to  study  that  still  need  to  be  considered.  One 
area  is  the  completeness  of  the  testing.  This  algorithm  was 
not  so  much  interested  in  computational  elegance  as  it  was 
in  correctness.  Investigating  the  long  run  times  and  trying 
to  cut  them  down  would  allow  more  complex  models  to  be  run 
without  the  extremely  long  waits .  One  thing  that  would 
certainly  help  would  be  to  put  this  algorithm  into  C  or 
Fortran  90 .  The  user  friendly  environment  of  MathCad 
doesn't  help  its  efficiency.  Also,  restructuring  the  sub¬ 
routines  and  cutting  down  on  redundant  calculations  will 
help  with  more  complex  models. 

Secondly,  the  formulation  of  the  model  itself  should  be 
considered  for  research.  The  ergodic  assumption  we  made  is 
not  true  for  a  branch  of  HMM  called  left  to  right  models. 
Once  the  state  is  transitioned  from,  it  is  never  allowed 
back.  Perhaps  this  type  of  model  would  better  suit  this 
application. 

Thirdly,  experimentally  validating  should  be 
considered.  Process  noise  and  measurement  noise  will  now 
play  a  role.  This  could  serve  to  ruin  the  model  altogether 
or  end  up  being  inconsequential.  The  algorithm  won't  serve 
any  purpose  if  it  doesn't  even  match  what  we  are  trying  to 
model . 
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Lastly,  this  algorithm  should  be  considered  with  other 
techniques  as  a  way  to  predict  movement.  For  instance,  the 
expected  Markov  state  could  be  used  as  a  way  to  measure 
noise  for  a  Kalman  filtering  system.  The  Kalman  filter 
would  then  predict  movement  based  on  the  constraints  imposed 
by  the  most  likely  Markov  state  [10] . 


CHAPTER  VII 


CONCLUSION 

Based  on  computer  generated  data,  we  used  HMM  to  try  to 
formulate  a  model  to  predict  future  movement  trends .  The 
useful  outputs  of  the  model  were  the  A,  B  and  %  matrices,  y, 
the  expected  state  at  time  t,  and  the  predictions  for  future 
movements.  All  outputs  are  based  on  bearing  only.  With 
velocities  and  time  intervals  for  specific  applications 
positions  can  be  inferred. 

MathCad  had  trouble  dealing  with  the  amount  of  small 
numbers  presented  to  it  so  numerical  stability  was  an  issue 
at  a  couple  of  different  places.  In  spite  of  that,  we  feel 
that  this  formulation  has  merit  and  should  be  considered 
when  implementing  movement  algorithms.  The  ability  to  model 
unseen  Markov  states  can  prove  to  be  very  powerful  in 
modeling  and  predicting  movements. 
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APPENDIX  A:  ESTIMATE  INITIAL  PARAMETERS  OF  MODEL 


This  model  is  for  4  possible  observation  symbols  and  3  possible  hidden  states.  Throughout 
the  text  this  model  will  be  referred  to  as  the  fundamental  model  because  it  is  the  basic  model 
that  we  will  use  to  expand  into  others. 


N;=14  i;=0..N  j:=0..N 


a(x,y)  :=.5-e  L 


(9x-8r+(9.y-4)^ 


These  functions  describe  each  of  the 
hills  in  the  z-function  below. 


b(x,y)  ;=.85-e 


’(9x-5)2+(9y- 8.5)2 
4 


General  form  of  these  equations  was  taken 
out  of  the  MathCad  5.0  reference  manual. 


f(x,y)  ;=.34-e  L 


'(9x-4.5)2+(9y-2)2 

■(9•x^-l)2  /'9•y^-l\l 

4  J 

-F.75*e 

1 

O 

+  b(x,y)-t-a(x,y) 


M.  .  :=10-f(Xi,yj) 


X  dwn  y  across 

In  the  above  contour  plot,  height  will  be  proportional  to  high  probabilities  in  the 
conceal  state,  where  as  the  inverse  will  be  proportional  to  high  probabilties  in 
speed  state. 
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In  order  to  start  at  gnd  zero,  we  need  to  difference  all  the 
elements  with  the  middle  one  where  we  start  the  vehicle  at. 


east  is  up, 
west  is  down, 
north  is  right, 
south  is  left. 


3.245 

2.836  2.453  2.092 

1.749  1.424 

1.118 

0.832 

3.04 

2.662 

2.308  1.966 

1.63  1.302 

0.99 

0.701 

2.782 

2.47 

2.178  1.874 

1.541  1.189  0.845 

0.536 

2.524  2.34 

2.17  1.935 

1.593  1.17 

0.74 

0.371 

2.312 

2.34 

2.369  2.242 

1.871  1.313 

0.723 

0.243 

2.138 

2.437 

2.723  2.734  2.323  1.585 

0.782 

0.156 

1.91 

2.458 

2.985  3.131 

2.698  1.813 

0.829 

0.08 

1.503 

2.179 

2.836  3.082  2.693  1.801 

0.781 

0 

4.282  4.92 

5.565  5.876 

5.647  4.964  4.116  3.411 

3.506 

3.993 

4.544  4.935 

5.005  4.743 

4.246  3.674 

2.744 

3.071 

3.544  4.077  4.558  4.84 

4.769  4.311 

2.109  2.335 

2.791  3.51 

4.409  5.209  5.519 

5.121 

1.625 

1.806 

2.27  3.12 

4.292  5.422 

5.968 

5.599 

1.256 

1.407 

1.84  2.676 

3.863  5.032 

5.628 

5.301 

0.961 

1.073 

1.412  2.081 

3.039  3.991 

4.484  4.228 

These  numbers 
correspond  to  the  heights 
of  the  above  z-function. 


This  is  the  matrix  that  will  be  used  in  order  to  find  the  highest  area  in  each  direction.  The 
routines  below  will  find  the  weighted  sum  across  the  respective  look  angles.  The  routines  are 
unique  for  each  direction. 


west : 


x<-0 

while  x^l4 
max^<-0. 

X4— x-h  1 
a<—4 
b^9 

while  a^O 

max^<-M^  ^  if  I  |>|  max^^ 
a«-a-  1 
max^«-inaxjj-.l 

a«“5 

b4-8 

while  a>0 

max,4-M^ ,  if  I  V  |>i  max, 
b  a,b  j  a,b  |  |  b 

a<-a-  1 

maXjj^-maXjj-.2 

a^— 6 
b<-7 

while  a>0 

max^,-M^  j,  if  I  |>|  max^ 
a<— a-  1 
max^<-max^j-.4 

a^-S 


b«— 6 

while  a^O 

maXb,-M^_b  I  Ma,b  |>|  maXb  [ 

af-a-  1 
maXj^4“  max^*.2 

a4-4 

b^5 

while  a^O 

max,4-M^  if  I  ,  |>|  max.  I 
b  a,b  j  a,b  |  |  b  | 

a<-a-  1 

maX|^4-max^j-.l 

14 

max<-  ^  maxj 
j  =  0 

max<-.5  if  max<0 
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east := 


x<— 0 

while  x^l4 


max^<-0 

X4-X+  1 

a^-lO 

b^9 

while  a<14 

max,+-M  .  if  I  M  ,  |>|  max, 
b  a,b  j  a,b  I  I  b 

a<-a-h  1 

maxjjf-max^j*.l 


a<-9 

b^8 

while  a<14 

maXjj«“M^  if  |  ^  j>j  maxj^ 

a<-a-|- 1 

maXjj<-maXjj-.2 

a<—  8 
b^7 

while  a^l4 

max.  ,  if  I  M  .  1>|  max, 

b  a,b  I  a,b  [  |  b 

a<-a-h  1 
max^4—  maXj^-.4 

a<-9 


b^6 

while  a^l4 

maXb,-M^_b  I  Ma,b  |>|  maXb 

a<-a-h  1 

maxjj+-max^*.2 

a^-lO 
b<— '  5 

while  a<14 

maXb^-M^  b  if  |  ^a.b  |>| 

a<-a-f- 1 

maxjj«~maxjj-.l 

14 

max<“  ^  maXj 

j  =  0 

max4-.5  if  max<0 
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south : 


jX4-0 

while  x^l4 

maXj^4-0 

x<-x-|- 1 
b4-4 
a<-9 

while  b>0 

maXa«-M^  b  I  ^a,b  |>| 

b*-b-  1 

max„4-max„'.l 
a  a 

b«-5 

8 

while  b>0 

maXa«-Ma  b  if  |  b  |>|  max^ 
b«-b-  1 
inax^<-max^-.2 

a  a 


a4-7 


while  b>0 

maXa^-Ma  b  if  I  Ma  b  |>|  maXa 

b«-b-  1 

niax^f-niax„-.4 
a  a 

b<-5 

a«-6 

while  b>0 

niaXa«-Ma_b  if  I  Ma  b  |>|  maXa 

b^-b-  1 

max„<“niax^-.2 
a  a 

b«-4 

a<-5 

while  b^O 

maXa«-Ma_b  if  I  Ma_b  |>|  maXa 

b«-b-  1 

max^^-max^-.l 

a  a 


14 


max4- 


max, 


j  =  0 

max4-.5  if  max<0 


north  := 


X4-0 

while  x<14 
jmax^<-0 


jx<-x-|“  1 
b4--10 
a^~9 


while  b^l4 

™“a-^a.b  if|Ma.bh™“a 

b<-b«i-l 

max„<-max„\l 

a  a 

b*“9 
a< —  8 

while  b^l4 

maXa^Ma  b  if  |  Ma_b  |>|  maXa 

b^—b  +  1 

max.«-max,  .2 

a  a 

b«-8 

a«-7 

while  b<14 

maXa«-Ma  b  if  I  Ma  b  |>|  maXa 
b«-b+l 
max„*“max„-.4 

a  a 

b*-9 
a<— 6 

while  b^l4 


"’“a^^a  b  if  I  Ma  b  |>|  maXa 

b*-b+l 

max^<-max^-.2 

a  a 

b^-lO 

a*-5 

while  b^l4 

™“a^^a,b  I  ^^a,b  |>|  ™“a 
b*-b+l 
max„*-max,-.l 

a  a 


14 


max<- 


maxj 


max<~.5  if  max<0 


west  =  0.844 


east  =  5.205 
north  =  7.757 
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south  =  3.808 

totconceal  := north  4-  east-f-  south -i-  west 


totconceal  =  17.614 


totconceal  will  be  used  as  the  normalizer  once 
all  the  respective  scores  are  added  up. 


for  conceal  state: 


north 

totconceal 


=  0.44 


west 

totconceal 


0.048 


A 


east 

totconceal 


0.296 


south 

totconceal 


=  0.216 


for  speed  state: 


totspd  := north'  ^  east'  ^  -f  south"  *  -|-  west" ' 


totspd  =  1 .769 


north"* 

totspd 


0.073 


west'  * 
totspd 


0.67 


A 


east'* 

totspd 


0.109 


south'  * 
totspd 


0.148 


SO,  we  can  see  how  the  "scores"  for  each  direction  add  up  to  form  probabilities  for  the 
respective  Markov  state. 
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let  0  be  the  conceal  state  and  2  be  the  speed  state.  North:0  east:1  south:2  west:3 


north  east 

totconceal  totconceal 


south  west 
totconceal  totconceal 


B  ;= 


.25 

north’  * 
totspd 


.25  .25  .25 

easf  *  south’  ’  west’  * 
totspd  totspd  totspd 


'0.44 

0.296 

0.216 

0.048  ■ 

B  = 

0.25 

0.25 

0.25 

0.25 

0.073 

0.109 

0.148 

0.67 

This  will  be  the  Initial  B  matrix  we  use  to  train 
the  model.  Notice  the  uniform  distribution  in 
the  ambiguous  state. 


w  := 0 ..  2  Define  range  variable  to  make  sure  the  rows  of  these  estimates  sum  to  unity. 


3 


j  =  0 


1 

T 

T 


summing  the  B  matrix  across  the  rows  yields  one.  This  iets  us  know  that 
the  matrix  is  vaiid  to  use. 


M  :=4  define  the  number  of  observations  possible. 

N  :=3  This  is  the  number  of  Markov  States  allowed, 
i  :=0..N-  1 


a  :=3 
p;=3 


Define  the  beta  distribution  and  the  parameters  we 
wiii  use. 


f(x,a,P) 


._  r(a  +  p) 

■  TT^OtTFT 


■(l-x)P- 


1 


j  ;=o..N-  1 
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Now  we  have  to  find  the  state  transition  matrix  (A)  first  guess.  Beta  distribution  wiil  be  used 
because  it  normalizes  to  one  and  is  easily  parameterized.  We  only  need  three  integrations 
because  there  are  only  three  states  to  transition  to. 


j.1 

N 

Jo 


f(x,a,P)dx  =  0.21 


f(x,a,P)dx  =0.21 


j; 


f(x,a,P)dx  =0.58 


Here  we  have  an  example  to  make  sure  we  are  doing 
things  properly. 


This  is  a  more  generalized  form  of  the  integrations  above.  Another  sample  calculation  tells  us 
that  we  are  still  okay. 


^(v,a,P) 


f(x,a,P)dx 


^(3, a, P)  =0.21 


■^(1,2,4)  ^(2,2,4)  ^(3,2,4)- 

’0.539  0.416  0.045' 

A:= 

^(1,3,3)  ^(2,3,3)  ^(3,3,3) 

A  = 

0.21  0.58  0.21 

_^(1,2,3)  ^(2,2,3)  ^(3,2,3)_ 

0.407  0.481  0.111 

A(a,b,c,d,e,f,N) 


^(N-2,a,b) 

4(N-2,c,d) 

,^(N-2,e,f) 


4(N-l,a,b)  ^(N,a,b) 
4(N-l,c,d)  ^(N,c,d) 
^(N-l,e,f)  ^(N,e,f) 


All  the  above  manipulation  was  either  putting  it  into 
an  easy  form  MathCad  will  agree  with  or  checking 
things  out  with  myself  that  they  actually  work. 
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A(2,4,3,3,4,2,3) 


0.539 

0.416 

0.045’ 

0.21 

0.58 

0.21 

0.045 

0.416 

0.539 

This  is  the  form  of  the  A  matrix  that  we  wiii  now  work  with  as  an  initiai  matrix 


A:=A(2,4,3,3,4,2,3) 


■  0.539 

0.416 

0.045  ‘ 

A  = 

0.21 

0.58 

0.21 

0.045 

0.416 

0.539 

notice  we  set  up  the  initiai  matrix  so  the  chance  of 
reentering  the  same  state  is  the  highest. 


j  =  0 
1 

T 
T 


W,J 


summing  the  A  matrix  across  the  rows  yieids  one.  This  iets  us  know  that 
the  matrix  is  vaiid  to  use. 


X  ;=0, .05..  1  define  a  range  variabie  to  piot  the  beta  functions  beiow. 


Plot  of  Beta  Functions  for  each  State 


here  we  have  the  initiai  beta  functions 
with  parameters  for  the  state 
transition  matrix.  These  were 
parameterized  this  way  because  they 
estimate  recurrent  states  being  more 
probabie  than  non-recurrent  states. 
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Now  the  initial  state  probability  matrix  {n)  has  to  be  determined.  It  is  the  same  as  the  state 
transition  matrix  but,  it  only  has  one  row. 


a  :=2 
p:=2 


f(x,a,P) 


r(a  +  p) 

-iTS[wm 


Show  the  beta  function  and 
parameters  again. 


j.1 

N 

Jo 


f(x,a,p)dx=0.259 


f(x,a,p)dx  =0.259 


1 

Jn 


f(x,a,p)dx  =0.481 


Again,  here  we  have  some  sample  calculations  to 
make  sure  the  equations  are  set  up  properly  and 
the  general  form  below  it. 


^(v,a,P) 


f(x,a,p)dx 


^(3, a, P)  =0.259 


7i;=(^(1,2,2)  ^(2,2,2)  ^(3,2,2)) 

7t=[  0.259  0.481  0.259] 


7t(a,b,N):=(4(N-2,a,b)  ^(N-l,a,b)  ^(N,a,b)) 


7t(2,2,3)  =[  0.259  0.481  0.259]  Here  we  set  the  most  probable  state  to 

start  in  as  the  ambiguous  state. 


n  ;=7t(2,2,3)  Again,  the  above  manipulations  are  just  to  get  the  notation  in 
a  form  that  is  easy  to  use  later  on. 


s  = 

j-0 


X  :=0,.05..  1 


2 


(f(x,2.2)) 


0 
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summing  the  n  matrix  across  the  row  yields  one.  This  lets  us  know  that 
the  matrix  is  valid  to  use. 


Plot  of  Beta  Function  for  Initial  State 


here  we  have  the  beta  function  with 
parameter  for  the  initial  state 
probability  matrix. 


APPENDIX  B:  Generate  Data 


The  T  tells  us  how  long  each  observation  set  is,  K  tells  us  how  many  observation  sets  to  iook 
at,  and  Reai  tells  us  how  many  sets  of  real  observations  we  have. 


T:=3  K:=2  Real:=5 

t:=0.. T  k:=0.. K  real  :=0.. Real 


time  :=0.. 2 


The  function  beiow  wili  give  us  random  numbers  but  they 
are  aii  based  on  the  beta  distribution.  The  parameters  are 
a  creative  choice. 


The  time  range  variable  is 
used  below  in  finding  the 
order  of  the  model. 


:=3.0rbeta(T-t-l,r,s) 


r:=2  s:=2 


We  are  only  interested  in  integers.  That  is  why  must  round  the  number  generator  to  make 
the  data  useful. 

round(x)  :=if(x-  floot(x)<.5,floot(x),ceil(x)) 


The  real  observation  data  wiii  be  used  later  on  to 
determine  the  order  of  the  modei.  The  obs  data 
taken  above  is  the  data  used  to  train  the  modei. 


0,k:=round(obs,^) 


realobs^“’^  :=3.5  rbeta(time-t-  l,r,s) 


re^lOtime.rear=rO“"‘l(’'ealobs^n,e,real) 


2  1 
1  2 
1  1 
1  1 


2‘ 

1 

3 

2 


realO  := 


‘1 

1 

2 


2 

0 

1 


2 

3 

2 


2 

1 

3 


2 

1 

2 


r 

3 

2 


realO  has  to  be  three  elements  long.  This 
is  true  because  it  is  what  we  use  beiow  to 
determine  the  order.  That  part  of  the 
algorithm  uses  3  eiements  exciusiveiy  to 
find  movement  trends. 


APPENDIX  C:  Forward  and  Backward  Variable 


FwdProb(col)  := 


i^O 
j^O 
set<— col 
while  1 

(o<set>)^ 

if-i-K  1 
while  1 

t4-0 

while  1 

N-  1 
suni4-  Y, 


1 

j*-j+  1 


s  =  0 
*sum*B 


.  frt<set>\ 


a 


Let's  look  at  just  the  last  column  of  data  from  Obs  matrix  to  show  this  variable.  As  time 
progresses  across  the  rows  we  can  see  how  small  the  probabilities  get  for  the  observed 
data. 


0.056 

0.017 

2.694*10"^ 

4.292*10'^ 

0.12 

0.027 

5.715*10"^ 

1.109*10"^ 

0.038 

5.273*  10”^ 

1.014*10“^ 

2.029*  10"^ 

FwdProb(O)  = 
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BackProb(col)  :=  i<-0 

set<-  col 
while  1 

Pj_T«- 1 

i«— i+1 
t*-T-  1 
while  t^O 
i«-0 

while  i^-  1 
N-  1 

j5) 

1 

tf-t-  1 

p 


Again,  let's  look  only  at  the  first  column  of  data.  The  sample  matrix  shows  us  starting  at 
one  and  working  back. 

■  0.017  0.067  0.268  1  ‘ 

BackProb(O)  =  0013  0.054  0.23  1 

8.724«10“^  0.038  0.176  1 
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As  time  increases  we  will  have  problems  with  underflow.  This  can  be  fixed  by  using  scaling 
factors  below.  Both  methods  should  produce  the  same  results. 


set  ;=0  specify  a  column  so  we  can  compare  apples  to  apples  below. 


.  these  two  functions  are  the  same,  they  just  are 

scale.  ;=-rT-T _ - _  set  up  differently. 


FwdProb(set)j  j 

i  =  0 


c(col)  := 


set*-  col 
t*-0 

while  t^T 


■  4.653  ■ 

■  4.653  ■ 

20.215 

c(0)  = 

20.215 

106.122 

106.122 

_  574.347  _ 

_  574.347  _ 

- 

y]  FwdProb(set)j  , 

t«-t+  1 


Two  methods  were  used  to  check  work.  Both  are  the  same  so  we  can  use  either  one  to  fix 
underflow  problems  later  on.  Now,  let's  add  them  to  the  Forward  and  Backward  variables. 


FwdProbhat(col)  l= 


set*-  col 
i*-0 


while  i^-  1 
t*-0 

while  t^T 

FwdProbhat;  ,*-c(set)j-FwdProb(set)j  ^ 
t«-t+  1 
i«-i  -t- 1 
FwdProbhat 


Hat  notation  indicates  that  the  scale  has  been  added  in. 


■  0.261 

0.342 

0.286 

0.247  ■ 

FwdProbhat(O)  = 

0.56 

0.552 

0.607 

0.637 

0.179 

0.107 

0.108 

0.117 

scale  term  restores  the  magnitude  of  the 
sum  downward  to  1 . 
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BackProbhat(col)  := 


set<-  col 

while  1 

t«-0 

while  t^T 

BackProbhatj  ^<-c(set)^*BackProb(set).  ^ 
1 

i<-i-l- 1 
BackProbhat 


Again  we  used  the  same  column  of  data.  These  numbers  go  well  over  one  because  the 
scale  is  based  on  the  forward  variable.  For  implementation  purposes,  it  is  better  to  use 
these  large  numbers  than  the  very  small  numbers  that  might  cause  underflow. 


■  0.077 

1.364 

28.456 

574.347  ■ 

BackProbhat(0)  = 

0.06 

1.091 

24.395 

574.347 

0.041 

0.764 

18.66 

574.347 

Now,  we  should  look  at  what  Markov  state  is  most  likely  at  any  particular  time. 
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Yi, 


^FwdProb(u)j  j-BackProb(u)j 
^FwdProb(u)j  j-BackProb{u)j 

i  =  0 


u  :=  0  u  is  the  column  of  data  in  the 
Obs  matrix  we  specify  for 
example  purposes. 


here  we  can  see  the  most  likely 
state  at  each  time.  Slightly 
different  optimization  criteria  than 
the  Viterbi  but  I  like  it. 


N-  1 


1 

T 

T 

T 


Yi.t 

’0.33 

0.406 

0.326 

0.247  ■ 

y  = 

0.551 

0.524 

0.593 

0.637 

0.119 

0.071 

0.081 

0.117 

only  one  of  the  columns  of  observations  was  picked  here  for  example 
purposes,  but  they  all  add  up  to  one  so  that  gives  us  a  good  feeling  that  all 
is  right  in  the  way  the  multiple  sequences  is  set  up. 


As  an  additional  check,  look  at  the  y  variable  again  with  the  scale  added  to  it.  It  should 
output  the  same  matrix. 


^i.t 


^FwdProbhat(u)j  j  BackProbhat(u)j 

■m - ^ ^ — 

y]  pwdProbhat(u)j  j  BackProbhat(u)j 
i  =  0 


N-  1 


i  =  0 


1 

T 

’0.33 

0.406 

0.326 

0.247  ■ 

y  = 

0.551 

0.524 

0.593 

0.637 

T 

T 

0.119 

0.071 

0.081 

0.117 

SO,  now  with  the  scaling 
constants  added  we  see  that 
they  cancel  out  each  other  and 
give  us  the  exact  answer  above. 
It  worksi 
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at  0: 


HatZero-  j 


BackProbhat(u)j  ,  (o<u>) 


•Fw<lProbhat(u)-  q 


N-  1  N-  1 


FwdProbhat(u).  Q*BackProbhat(u)j  j*Aj  j-B  / 
i^o  j  =  0  '  '  ’ 


N-  1 


s 

j  =  0 
0.33 

ir53T 

UTT? 


HatZerOj  j 


^i.o 

0.33 

I533T 

irm 


t 


T 


c,  :=  P]  scale^ 
T  =  0 


:=  pj  scale^ 

T  ==  t 


We  will  use  a  different  way  than  matrix  algorithm  to  check  the  output. 


TestA.  .  := 


T-  1 

S  (o<o>)^^j  BackProb(0)j 


T-  1  N-  1 


.S  (Q<o>)^^j-BackProb(0)j 


TestA  = 


'0.629  0.358 
0.303  0.618 
0.094  0.627 


0.014' 

0.079 

0.279 


after  the  matrix  routine  calculates  the  two  A  matrices 
should  be  the  same.  They  are  in  this  case  to  about  3 
one-thousandths.  That  is  close  enough  for 
Government  work.  If  we  were  to  shorten  the  time 
interval  by  one  the  elements  would  be  totally  equal. 
MathCad  begins  to  fail  us  numerically  with  bigger  time 
step. 


This  is  part  of  the  parameter  re-estimation.  A  new  variable  has  to  be  defined  to  make  sure 
we  don't  go  out  of  bounds  in  the  re-estimation  algorithm. 


T  ;=0..T-  1 
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Matrix  : 


K  := 0  Here  we  use  only  the  first  column  because  we  want  to  check  by  doing  it  a 
different  way  with  the  TestA  above. 

u<-0 

while  u^K 
T^O 
i^O 
j^O 
x«-0 
y^O 

while  x^-  1 
while  y:^-  1 
sumZeta^  ^<-0 

y«-y+l 
y^O 

X<-X-|-l 

while  1 

while  i^-  1 
while  j  ^  “  1 

BackProbhat(u)j  j-B^  ^-FwdProbhatCu); 

“i.j'  N-  1  bJ-  1 

Y  FwdProbhat(u)  ■BackProbhat(u)b  ,-A  J,-B  .  . 

a=0  b  =  0 

sumZeta.  sumZeta.  .  +•  H-  . 


i.j 


ij 


j«-j+ 1 

i<-i  -h  1 
j^O 
T4— T  1 
i^-0 

while  t<T 

while  i^-  1 

^FwdProbhat(u).  ^•BackProbhat(u)j 

: 


Yi.t" 


^FwdProbhat(u).  ^‘BackProbhat(u). 

i  =  0 


i<-i-h  1 

t4-t+  1 

i<-0 
while  i^-  1 


suniGaniina.+ 


T-  1 


V  Y. 


Yi.t 


sumGammaj 

1 

-0 


while  x^“  1 
while  y^"  1 
Abar  <-0 

y«-y+  1 
y«-0 
x<— X-f  1 
while  1 

I  while  i^-  1 


sumZeta.  ■ 

Abar.  .< - ^ - — 

*»J  sumGamma. 


jH  + 1 

i<-0 

i«-0 

while  1 

i^i  ^  I 

i<-0 

while  i^-  1 

newsumGamma.f 

i<“i-f  1 
t«-0 


T 

Sr 

t-0 


y<-0 

while  x^-  1 
[while  y<M"-  1 


[y^y-bi 

y<-0 
X«“X*f  1 

while  i^-  1 
I  while  t^T 
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i.  o  , 
1 

i«— i  -I- 1 
t<-0 
i«— 0 
j^O 
£  <—  .02 

while  1 

while  j  ~  1 


i,  o 


Bbar. 


Yij 


newsumGamma. 


Bbarj  jf- 

jH-b  1 


'8  if  Bbar-  j<6 


M-  1 


normalizer-f- 

j^O 

while  j  - 


Bbar. 


L-0 


Bbar. 


i,L 


Bbar. 


i.J 


normalize 


j<-j  -I- 1 
i+l 
0 

B<—  Bbar 
7t4-7rbar 
A^-  Abar 
u-f  1 

x<-0 

y<-0 

while  x^-  1 

while  y<2-N+M-  1 


out„ 


-0 


y^-y-b  1 

y<-0 
Xf-x-h  1 
x<“0 
y<-0 

while  x^“  1 
while  1 

y«-y-hl 
y*-0 
x<— X*f  1 
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x<-0 

Vf-N 

while  x^-  1 

while  1 

0“tx,y«-Bx.(y_N) 

y*-y+i 

y<-N 
X<-X-|-  1 
x<-0 

y«-N+-M 

while  y:^N-HM-|-N“  1 
°'*^x,y*~’^0,(y-(N-|-M)) 

y*-y+ 1 

out 
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Here  we  output  all  the  matrices  as  an  aggregate  because  this  is  the  way  MathCad 
programming  structures  work.  The  programming  structures  below  it  simply  extract  the  proper 
matrix  form  the  aggregate. 


10.626  0.36  0.014  0.019  0.719  0.242  0.019  0.33  0.551 

03  0.619  0.081  0.019  0.732  0.23  0.019  0  0 

0.093  0.625  0.282  0.019  0.665  0.296  0.019  0  0 


A  := 


i<-0 

while  1 


A^' Matrix"^' ^ 


if-i-h  1 


B  := 


while  1 

1 

B 


■  0.626 

0.36 

0.014' 

■  0.019 

0.719 

0.242 

0.019' 

A  = 

0.3 

0.619 

0.081 

B  = 

0.019 

0.732 

0.23 

0;019 

0.093 

0.625 

0.282 

0.019 

0.665 

0.296 

0.019 

2 


1 

T 

1 


Rows  sum  to  one  here! 


3 


1 

T 

1 


Rows  sum  to  one  here 
also! 


71  := 


i4„N+.M 

while  i^N-t-M-i-N-  1 
i4-i-|-l 


7C 


7r=[0.33  0.551  0.119] 

2 

^  tcq  j  =  1  The  row  sums  to  one  here! 
j  =  0 


Since  all  the  rows  sum  to  one  for  a  the  matrices,  we  can  say  that  these  are  valid  to  work  with 
now. 


APPENDIX  D:  Most  Probable  Sequences 


We  need  to  redefine  the  FwdProb  and  scale  here  because  we  are  working  with  new  matrices. 

The  previous  routines  worked  with  the  estimates  and  these  will  work  with  the  final  output  matrices. 


col :=  2 


FwdProb(col)  := 


i«-0 

j«-0  c(col) 

set<—col 

while  1 

(o<set>)^ 

1 

while  1 

t<-0 

while  t^T-  1 

N-  1 

sum«-  >  a,  • 
s  =  0 

t«-t-|-l 

jH  + 1 

a 


set*-  col 
t*-0 

while  t^T 


"“•m - 

FwdProb(set)j  ^ 

i^ 

t*-t+l 


n  :=4  possible  observations  u  •  ^  r  i.  • 

^  changing  the  number  of  observations  is  easy, 

but  the  length  will  require  some  work.  We 

._  3  would  have  to  add  another  nested  loop  to  an 

length  of  observations  already  hairy  code. 

The  sequence  algorithm  enumerates  all  the  possible  sequences.  We  are  interested  here  in  only 
three  time  steps  so  we  will  redefine  the  maximum  time  value  at  2.  Remember  that  MathCad 
starts  its  elements  at  0  instead  of  1 . 


sequence : 


sequence 


q«-“0 

i^-O 

while  i<r 

Pi.o^O 

if-i-i- 1 
if-0 


while  j<n^ 
while  i<r 
sea  .4-0 

i4-i-f  1 
i<-0 
jH-hi 
i<-0 
j^O 

while  Pj._3<n 


while 

while  Pj._  j  <n 
seq^^<-p 


|q<-q+i 

Pr-l^O 

Pr-2*-Pr-2+l 


Pr-2^0 

Pr_3^Pr-3+l 


Pr-3^Pr-3“  ^ 

q 


T:=2 


64  n'^  =  64 


sequence  equals  what  combinatorics  predicts  so  this  routine  is  correct. 


LOOP  routine  calculates  the  probability  of  each  sequence.  It  will  be  used  as  sort  of  a 
subroutine  in  the  next  structure. 


LOGP(p) : 


setf-0 

i«-0 

j<-0 

while  1 

1 

while  j  ^  -  1 
t<-0 

while  t^T-  1 
N-  1 

sum^  X!  “s.tA.j 

s=  0 

1 

j<-j  -I- 1 

N-  1 

probab<—  ^  Oj  ^ 


probab 


Bestset  := 
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q^—O 
j^O 
norm<— 0 


numcol*- 

i4-0 

while  i<numcol 
totj<-0 
1 

if-0 

while  i<r 


i<-i^  1 

while  j<numcol 
while  i<r 
max-  j<— 0 

i«-i-h  1 
i^O 
jH-f  1 
i<-0 
while 


while  Pj._2<n 
while  p^_j<n 

probability^  LOGP(p) 

max'^^^-p 

totq<~  probability 

Pr-l^Pr-l+l 

q«-q+ 1 

Pr-l"“® 

Pr-2^Pr-2+l 

Pr-2*”® 

Pr-3^Pr-3+l 
Pr-3^Pr-3-l 
answer^  Q<-q 


numcol-  1 


normal*- 


E 

i  =  0 


answer^  j<- normal 
i*-0 

while  i<numcol 


bestset  uses 
the  logP 
routine  to  find 
each 

sequences 
probability  and 
then  rank 
order  them 
appropriately 


while  i<numcol 


totj 

answer,  - r 

^ » *  normal 

1^2 

j^O 

while  j<numcol 
while  i<r-h2 

answer^  j«-max. 

i<— i*h  1 
i<-2 

jH  +  i 

answer 


bestset  :=rsort(Bestset,  1 ) 
T 

bestset  := bestset 

bestset : = reverse(  bestset ) 
T 

bestset  := bestset 


These  statements  sort  the  above  matrix  according 
to  the  probabilities.  Then  the  matrix  is  output  with 
the  highest  probability  first  etc. 


This  is  the  stacked  Bestset  output.  Here  we  can  see  what  movements  will  most  likely  occur 
over  the  next  3  time  iterations. 


Log  Probability 


APPENDIX  E:  Order  Calculations 


realO  = 


1 

1 

2 


22221' 
0  3  113 
1  2  3  2  2 


This  data  was  generated  above  in  Appendix  B.  It  will  now 
be  used  to  "score"  the  model.  We  can  find  the  probabilities 
of  these  data  sequences  using  the  LOOP  routine.  The 
routine  uses  the  final  output  matrices  from  the  training 
data.  If  the  matrices  are  perfect  then  the  probabilities 
should  be  close  to  zero  on  the  log  scale.  Therefore,  if  we 
add  up  all  the  log  probabilities,  the  lowest  score  will  be  the 
best  model.  This  will  be  because  the  matrices  were  closest 
to  being  perfect  for  that  formulation. 


sightings^,  :=log 


/LOGP(realO"^'>)\ 
1  Bestsetg  j 


sightings  = 


■-0.903' 

-2.476 

-2.954 

-2.475 

-1.379 

-2.478 


i  :=0.. Real 


Here  we  will  sort  the  sightings  from  highest  probability  to  smallest  probability, 
sightings  :=sort( sightings) 


sightings  :=reverse( sightings) 


sightings  = 


-0.903 

-1.379 

-2.475 

-2.476 

-2.478 

-2.954 


0 


1 


4 


5 
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Real 

Score  :=  ^  sightings j 
i  =  0 

Score  =-12.666 


Real  ;=63 


Real 

prob2  :=  ^  bestsetj  j 
i  =  0 


The  Rating  closest  to  zero  will  be  the 
best  model. 


For  the  64  possible  sequences,  the  probability  of  each 
event  summed  should  add  up  to  one.  After  normalization  it 
did. 


prob2  =  1 


APPENDIX  F:  Alterations  for  6  Observation  Symbol  Model. 

These  routines  are  used  for  each  direction  to  look  at  the  weighted  sum  of  the  highest  points 
over  that  particular  swath.  They  are  unique  for  each  direction. 


northwest  :=  x«- 1 

while  x^3 

niaXj^<-0 

X4-X+1 

a<-5 
b4-8 
c<—  I 

while  a^O 

inax^«-M^  u  if  I  K  |>|  max^  I 

c  a,D  I  a,b  I  I  c  I 
af-a-  1 

max  <— max^\5 

c  c 

a«~4 

b^8 

c«-2 

while  a^O 

max^4-M, .  if  I  u  |>|  max^  I 

c  a^D  I  a,D  ]  I  c  I 

a4— a-  1 

bf-b-hl 

max  ♦-max^-.lS 
c  c 

a«— 3 
b«-8 
c«-3 

while  a^O 

max„«—  M,  u  if  I  M.  1 1>|  max,  I 
c  a,b  I  a,b  |  j  c  j 

a«-a-  1 

b<-b+l 

max^<-max^*.25 

c  c 

3 

max4—  ^  maXj 

j  =  0 

max4-  .5  if  max<0 


southwest  := 


X4-  1 

while  x<3 
max^4-0 

x<-x-h  1 

a4-5 

b«—6 

1 

while  a^O 

max^«~M,  if  I  K  l>l 

c  a,b  I  a,b  I  I  c 
a*-  a  -  1 
b^b-1 
max  <— max  -.5 

c  c 

a*- 4 
b<-6 
2 

while  a^O 

max^<— u  if  1  M,  1,  |>|  max^ 

c  a,b  {  a,b  I  I  c 

a4-a-  1 

b^-b-  1 
max^<-max^*.25 

c  c 

a<-3 

b*-6 

c<-3 

while  a^O 

max^4“M^ .  if  I  V  |>|  max^ 

c  a,b  I  a,b  |  |  c 

a<-a-  1 

b<-b-  1 

max^4-max^*.25 
c  c 


3 


max«-.5  if  max<0 
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southeast : 


1 

while  x:^3 
max^4-0 

x<-x-t- 1 
a<~9 
b<-6 
c<-l 

while  a^l4 

max^«-M^_b  if  |  i^a,b  |>| 

a4-a-}- 1 
bi-b-1 
max  <-“max  -.5 

c  c 

a^lO 

b4-6 

c^2 

while  a^l4 

max^<-M„  u  if  I  .  |>|  max^ 

c  a,b  I  a,b  |  |  c 

a<— a-h  1 
b<-b-  1 
max^<— max  -.25 

c  c 

a^-ll 
b4— 6 

c<-3 

while  a<14 

max^<-M^  .  if  I  M„  u  l>l  niax 

c  a,b  I  a,b  |  |  c 

a<— a-h  1 
b<-b-  1 
max^«“max^-.25 

c  c 

3 

max<-  ^  maxj 
j  =  0 

max^.5  if  max<0 


northeast  := 


x<-i 

while  x^3 

maXj^<-0 

x<— x-h  1 


a<-9 


b4--8 
c<- 1 

while  a^l4 

max^f-M^  b  if  |  i^a,b  |>| 

a<— a-h  1 
b^b-hl 
max^<-max^*.5 

c  c 

a^lO 

b4-8 


C+-*  2 

while  a^l4 

max^«-M^  b  if  I  Ma,b  |>1  max^ 

a*-a-f- 1 
b<-b-|-  1 
max^<-max^'.25 

c  c 

a<-ll 
b<-8 
Cf-3 


while  a<14 


max^ 


•  .  if  M, ,  >  max^  | 

a,b  a,b  c 


a*“a-f  1 
b4-b-h  1 


max^«“max^*.25 

c  c 


maXj 


max4—  ^ 

j  =  0 

max4-.5  if  max<0 


ssouth  : 


X4-  1 

while  x!^3 

maXj^4-0 

x<— x-h  1 
2i^l 

b<-6 
c<- 1 

while  b^O 

max^^M^  b  I  Ma,b  |>i 

b4“b-  1 

max^^max^*.5 
c  c 

a<“6 

b<-6 

c<“  2 

while  b^O 

max^-M^_b  if  |  i>| 

b«-b-  1 
max  <-“max  -.25 

c  c 

a<-“8 
h^6 
c<-3 

while  b>0 

max^<-M^_b  if  i  Ma,b  |>|  max^ 

b<-b-  1 
max^4“max  -.25 

c  c 

3 

max«—  ^  maxj 
j  =  0 

max<— .5  if  max<0 


nnorth  := 
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x<~  1 

while  x<3 

Imax^<-0 

X4-X-I-  1 
a4^7 
b^8 

C4-  1 

while  b<14 

|max^<-M^_b  if  |  ^a.b  |>| 

jb-b-t-l 

max„f-max^-.5 
c  c 

a<— 6 

bf-8 

c<“  2 

while  b^l4 

Imax„4-M„  .  if  I  M  .  |>|  max 

C  .3,0  I  a,D  I  I  c 

b*-b-h  1 
max„<— max  -.25 

c  c 

a«-8 

b<-8 

c<-3 

while  b<14 

max^-M^  b  if  |  ^a.b  |>| 

b«-b-|- 1 
max^4-max  -.25 

c  c 

3 

maxf-  ^  maXj 
j  =  0 

max<-.5  if  max<0 


northwest  =0.5 
ssouth  =  3.792 
southeast  =  5.102 


southwest  =  2.554 
nnorth  =  8.339 
northeast  =  3.62 
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let  0  be  the  conceal  state  and  2  be  the  speed  state.  NNorth:0  northeast:  1  southeast:2 
SSouth:3  southwest:4  northwest:5.  Keep  In  mind  here  that  the  number  of  states  doesn't 
change  from  3.  That  is  why  there  are  3  rows  in  the  B  matrix  still. 


nnorth 

northeast 

southeast 

ssouth 

southwest 

northwest 

totconceal 

totconceal 

totconceal 

totconceal 

totconceal 

totconceal 

.166667 

.166667 

.166667 

.166667 

.166667 

.166667 

nnorth"  ^ 

northeast"  ^ 

southeast"* 

ssouth"* 

southwest"  * 

northwest"  * 

totspd 

totspd 

totspd 

totspd 

totspd 

totspd 

■  0.349 

0.151 

0.213 

0.159 

0.107 

0.021 

0.167 

0.167 

0.167 

0.167 

0.167 

0.167 

0.037 

0.085 

0.06 

0.081 

0.121 

0.616 

This  will  be  the  initial  B  matrix  we  use  to  train 
the  model.  Notice  the  uniform  distribution  in 
the  ambiguous  state. 


w  :=  0..  2  Define  a  range  variabie  to  sum  rows  of  the  B  matrix. 


5 


j  =  0 


1 

T 

T 


All  the  rows  still  sum  to  one. 


M  ;= 6  define  the  number  of  observations  possible. 


From  here  on  down  Appendix  A  wiii  pick  up  nicely  since  nothing  else  changes. 


APPENDIX  G:  Alterations  for  4  Markov  State  Model. 


The  number  of  symbols  reverts  back 
to  4  here  for  symplicity. 


f(x,a,P) 


r(a  +  p)  „_i 


•d-x) 


p-i 


N:=4 


The  beta  distribution  is  here  because  that  is 
how  we  plan  to  estimate  the  changed  states. 


4(v.a,P) 


(v-1) 


1 

'TT 


f(x,a,P)dx 


now  that  we  have  more  than  3  states, 
defining  will  be  allittle  more  difficult,  so, 
let's  just  use  the  beta  to  fill  in  the  blanks. 


let  0  be  the  conceal  state  and  2  be  the  speed  state.  NorthiO  east;1  south:2  west:3 


north 

east 

south 

west 

totconceal 

totconceal 

totconceal 

totconceal 

4(1.2.4) 

^(2,2,4) 

^(3.2,4) 

^(4.2,4) 

^(1,4,2) 

^(2,4,2) 

^(3,4,2) 

^(4,4,2) 

north"  ’ 

east"’ 

south"’ 

west"  ’ 

totspd 

totspd 

totspd 

totspd 

We  have  replaced  the  ambivalent 
state  with  2  semi-ambivalent 
states.  There  is  really  no  physical 
significance  that  we  can 
determine  now  but  that  doesn't 
matter.  The  parameters  are  set 
up  to  compliment  the  other 
distribution.. 


0.44 

0.296 

0.216 

0.048 

0.367 

0.445 

0.172 

0.016 

0.016 

0.172 

0.445 

0.367 

0.073 

0.109 

0.148 

0.67 

B.. 


W.J 


w:=0..N-  1 
All  rows  sum  to  one. 


This  will  be  the  initial  B  matrix  we  use  to  train 
the  model.  Notice  the  uniform  distribution  in 
the  ambiguous  state. 


M  ■.=4 


define  the  number  of  observations  possible. 


104 


N:=4 

i  :=0..N-  1 
j  :=0..N-  1 


a  :=3 
p:=4 


V  /V  ft  \  '^P)  Cl—l^l  

f[x,a,p)  ■"  ■*  '(i"*) 


Now  we  have  to  find  the  state  transition  matrix  first  guess.  Beta  distribution  will  be  used 
because  it  normalizes  to  one  and  is  easily  parameterized.  We  only  need  three  integrations 
because  there  are  only  three  states  to  transition  to. 


•  1 


Jo 


f(x,a,p)dx=0.169 


f(x,a,P)dx  =0.306 


1 

TT 


T7 

1 

T7 


fl:x,a,P)dx=0.487 


P[x,a,P)dx  =0.038 


^  ll 

V 


^(v,a,P):= 


TT 


(v-1) 


1 

'TT  . 


f(x,a,P)dx 


When  doing  sample  calculations  now 
an  extra  integartion  is  needed 
because  we  have  an  extra  transition 
possible. 


^(3, a, P)  =0.306 


A;= 


4(1, 2, 4)  4(2, 2, 4)  4(3,2,4)  ^(4,2,4) 
^(1,4,3)  4(2,4, 3)  ^(3,4,3)  ^(4,4,3) 
4(1,3,4)  4(2,3.4)  ^(3,3,4)  ^(4,3,4) 
^(1,4,2)  ^(2,4,2)  40,4,2)  4(4,4,2) 


A(a,b,c,d,e,f,g,h,N)  := 


4(N- 


■  0.367 

0.445 

0.172 

0.016' 

0.038 

0.306 

0.487 

0.169. 

A  = 

0.169 

0.487 

0.306 

0.038 

- 

0.016 

0.172 

0.445 

0.367 

4(N-2,a,b) 

4(N- 

l,a,b) 

4(N,a, 

,b)l 

4(N-2,c,d) 

^(N- 

l,c,d) 

4(N,c,d) 

4(N- 

2,e,f) 

4(N- 

l,e,f) 

4(N,e.f) 

4(N-2,g,h) 

4(N- 

l,g.h) 

^(N,g,h) 
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A(2, 4, 4, 3, 3, 4, 4, 2,4)  = 


0.367 

0.445 

0.172 

0.016" 

0.038 

0.306 

0.487 

0.169 

0.169 

0.487 

0.306 

0.038 

0.016 

0.172 

0.445 

0.367 

A:=A(2,4,4,3,3,4,4,2,4) 


All  the  above  manipulation  was  either  putting  it  into 
an  easy  form  MathCad  wil  agree  with  or  checking 
things  out  with  myself  that  they  actually  work. 


■  0.367 

0.445 

0.172 

0.016 

0.038 

0.306 

0.487 

0.169 

0.169 

0.487 

0.306 

0.038 

0.016 

0.172 

0.445 

0.367 

notice  we  set  up  the  initial  matrix  so  the  chance  of 
reentering  the  same  state  is  the  highest. 


3 


1 

T 

T 

T 


All  rows  still  sum  to  one  so  we 
know  this  is  a  valid  matrix. 


X  :=0,.05..  1 

Plot  of  Beta  Function  for  each  State 


here  we  have  the  initial  beta  functions 
with  parameters  for  the  state 
transition  matrix.  Now  there  is  simply 
an  extra  state  to  deal  with. 


106 

Now  the  initial  state  probability  matrix  has  to  be  determined.  It  is  the  same  as  the  state 
transition  matrix  but,  it  only  has  one  row. 


The  procedure  is  the  same  here 
as  in  Appendix  A,  basically.  The 
difference  is  the  extra  element  to 
account  for  the  extra  Markov  state. 


a  ;=2 


p:=2 

f(x,a,P) 


._  rfa  +  P) 
-TTSTTdiV^ 


•(l-x) 


p-i 


|.i 

N 

Jo 


f(x,a,p)dx  =0.156 


f(x,a,p)dx  =0.344 


^(v,a,P) 


f(x,a,P)dx 


1 

Jn 


fl[x,a,P)dx  =0.344 


f(x,a,p)dx  =0.156 


4(3,a,p)  =0.344 


7t  :=(4(1,2,2)  ^(2,2,2)  ^(3,2,2)  ^(4,2,2))  7t  =[0.156  0.344  0.344  0.156] 
7i(a,b,N):=(^(N-3,a,b)  ^(N-2,a,b)  ^(N-l,a,b)  ^(N,a,b)) 


71(2,2,4)  =[0.156  0.344  0.344  0.156] 


X  ;=0,.05..  1 


t:  :=:t(2,2,4) 


Plot  of  Beta  Function  for  Initial  State 


here  we  have  the  beta  function  with 
parameter  for  the  initial  state 
probability  matrix. 


APPENDIX  H:  Alterations  for  5  Markov  State  Model. 


rr  a  4-  R  number  of  symbols  is  going  to  be 

f(x,a,p)  := '•(  1  -  x)P  - '  4  here  also  for  simplicity  sake. 


N  :=4  here  N  is  just  the  number  of  integrations  we  have  to  do.  For  right  now,  it 
doesn't  stand  for  the  number  of  Markov  states  possible. 


^(v,a,P) 


f(x,a,P)dx 


now  that  we  have  more  than  3  states, 
defining  will  be  allittle  more  difficult,  so, 
let's  just  use  the  handy  dandy  beta  to  fill  in 
the  blanks. 


let  0  be  the  conceal  state  and  2  be  the  speed  state.  North:0  east:1  south;2  west:3 


north 

east 

south 

west 

totconceal 

4(1,2, 4) 
.25 

4(1,4, 2) 
north*  ^ 

totconceal 

4(2, 2,4) 
.25 

4(2, 4, 2) 
east*  * 

totconceal 

40,2,4) 

.25 

4(3,4, 2) 
south* ' 

totconceal 

4(4,2,4) 

.25 

4(4,4, 2) 
west*  * 

Here  we  added  a  true 
ambivalent  state  back  into  the 
model.  How  we  choose  to 
estimate  these  states  is  a 
creative  choice  and  this  is  what 
we  have  choosen. 

totspd 

totspd 

totspd 

totspd 

0.44 

0.296 

0.216 

0.048  ■ 

0.367 

0.445 

0.172 

0.016 

0.25 

0.25 

0.25 

0.25 

0.016 

0.172 

0.445 

0.367 

0.073 

0.109 

0.148 

0.67 

3 


w:=0..N 


E 

j-O 

1 

T 

T 

T 

1 


B 


W.J 


All  rows  still  sum  to  one. 


This  will  be  the  initial  B  matrix  we  use  to  train 
the  model.  Notice  the  uniform  distribution  in 
the  ambivalent  state. 


M  :=4  define  the  number  of  observations  possible. 


Now  we  let  N  be  the  states  allowed  in  the  model. 
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N:=5 

i  :=0..N-  1 
j:=0..N-l 


a  :=3 


p;=4 

f(x,a,P) 


r(a  +  P) 

■■r(a)r(p)-^ 


(l-x) 


p-i 


Now  we  have  to  find  the  state  transition  matrix  first  guess.  Beta  distribution  will  be  used 
because  it  normalizes  to  one  and  is  easily  parameterized.  We  only  need  three  integrations 
because  there  are  only  three  states  to  transition  to. 


•  1 

TJ 

Jo 


f(x,a,p)dx  =0.099 


f(x,a,P)dx  =0.365 


5- 


4 


f(x,a,P)dx  =0.017 


1 

N 


f(x,a,p)dx=0.357 


f(x,a,p)dx=0.162 


Notice  again  that  the  number  of  sample  calculations 
has  to  increase  by  one  because  of  the  extra  state. 


^(v,a,P) 


(V- 


f(x,a,P)dx 


^(3, a, P)  =0.365 


A:= 


b 


^(1,2,4)  ^(2,2,4) 

^(1,4,3)  ^(2,4,3) 

^(1,3.5, 3.5)  ^(2, 3.5, 3.5) 
^(1,3,4)  ^(2,3,4) 

^(1,4,2)  ^(2,4,2) 


^(3,2,4) 
^(3,4,3) 
^(3, 3.5, 3.5) 
^(3,3,4) 
^(3,4,2) 


^(4,2,4)  4(5,2,4)  ■ 

^(4,4,3)  ^(5,4,3) 

^(4,3.5,3.5)  ^(5, 3.5, 3.5) 
^(4,3,4)  4(5, 3, 4) 

4(4,4, 2)  4(5, 4, 2) 


0.263 

0.4 

0.25 

0.08 

6.72*  10' 

0.017 

0.162 

0.365 

0.357 

0.099 

0.044 

0.259 

0.394 

0.259 

0.044 

0.099 

0.357 

0.365 

0.162 

0.017 

6.72»10"^ 

0.08 

0.25 

0.4 

0.263 

A(a,b,c,d,e,f,g,h,N)  := 


■^(N-3,a,b)  ^(N-2,a,b) 
^(N-3,c,d)  ^(N-2,c,d) 
^(N-3,e,f)  ^(N-2,e,f) 
_^(N-3,g,h)  ^(N-2,g,h) 


^(N-l,a,b)  4(N,a,b)- 
4(N-l,c,d)  ^(N,c,d) 
^(N-l,e,f)  ^(N,e,f) 
^(N-l,g,h)  ^(N,g,h) 


A(a,b,c,d,e,f,g,h,i,j,N) 


■^(N-4,a,b)  ^(N-3,a,b) 
^(N-4,c,d)  ^(N-3,c,d) 
^(N-4,e,f)  ^(N-3,e,f) 
4(N-4,g,h)  4(N-3,g,h) 
J(N-4,i,j)  ^(N-3,i,j) 


^(N-2,a,b)  ^(N-l,a,b)  ^(N,a,b) 
^(N-2,c,d)  ^(N-l,c,d)  ^(N,c,d) 
^(N-2,e,f)  ^(N-l,e,f)  ^(N,e,f) 
^(N-2,g,h)  ^(N-l,g,h)  ^(N,g,h) 


A(2, 4, 4, 3, 3.5, 3.5,3, 4, 4, 2, 5) 


0.263 

0.4 

0.25 

0.08 

6.72*1 

0.017 

0.162 

0.365 

0.357 

0.099 

0.044 

0.259 

0.394 

0.259 

0.044 

0.099 

0.357 

0.365 

0.162 

0.017 

6.72*10“^ 

0.08 

0.25 

0.4 

0.263 

All  the  above  manipulation  was  either  putting  it  into 
A  ;=A(2,4,4,3,3.5,3.5,3,4,4,2,5)  an  easy  form  MathCad  wil  agree  with  or  checking 

things  out  with  myself  that  they  actually  work. 


0.263 

0.4 

0.25 

0.08 

6.72*1 

0.017 

0.162 

0.365 

0.357 

0.099 

0.044 

0.259 

0.394 

0.259 

0.044 

0.099 

0.357 

0.365 

0.162 

0.017 

6.72*10“^ 

0.08 

0.25 

0.4 

0.263 

notice  we  set  up  the  initial  matrix  so  the  chance  of 
reentering  the  same  state  is  the  highest. 
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4 

Z  ^v.j 

j  =  0 

T1 


All  rows  still  sum  to  one  indicating  that  the  matrix  is  valid. 


X  :=0,.05..  1 


fl:x,2,4) 
fix, 4, 3) 
■f(x,"3,4) 
'?x'[4,2) 
11x73.5,3.5) 


Plot  of  Beta  Function  for  each  State 


here  we  have  the  initial  beta  functions 
with  parameters  for  the  state 
transition  matrix. 


Ill 

Now  the  initial  state  probability  matrix  has  to  be  determined.  It  is  the  same  as  the  state 
transition  matrix  but,  it  only  has  one  row. 


The  procedure  is  the  same  here 
as  in  Appendix  A,  basically.  The 
difference  is  the  extra  two  elements  to 
account  for  the  extra  Markov  states. 


a  :=2 


p;=2 


f(x,a,P) 


r(a-t-P)  ^a-l 
■“r(a)r(p)-^ 


•(l-x)P-‘ 


•  1 

N 

Jo 


f(x,a,P)dx=0.104 


2 


f(x,a,P)dx  =0.296 


5- 


f(x,a,p)dx  =0.104 


f(x,a,P)dx=0.248 


fl[x,a,p)dx  =0.248 


^(v,a,p) 


(V-  1) 


lfx,a,p)dx 


^(3, a, p)  =0.296 


71  :=(^(1,2,2)  4(2, 2, 2)  4(3, 2, 2)  4(4, 2, 2)  4(5, 2, 2)) 
71  =[  0.104  0.248  0.296  0.248  0.104] 


7t(a,b,N);=(4(N-4,a,b)  4(N-3,a,b)  4(N-2,a,b)  4(N-l,a,b)  4(N,a,b)) 
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71(2,2,5) 
X  ;=0,.05..1 

(fl:x,2,2)) 


4 

E  ’^o.j 

j  =  0 


=  [0.104  0.248  0.296  0.248  0.104] 


n  :=n(2,2,5) 


Plot  of  Beta  Function  for  Intial  State 


here  we  have  the  beta  function  with 
parameter  for  the  initial  state 
probability  matrix. 


1 


All  rows  still  sum  to  one. 


APPENDIX  1:  Correctness  of  Model. 


Matrices  for  trivial  sequence  0,0,0 

0.72  0.273  7.346»10"^  T  0.943  0.019  0.019  0.019" 

^  ~  0.409  0.544  0.047  ®  “  0.943  0.019  0.019  0.019 

0.153  0.653  0.193  0.943  0.019  0.019  0.019 


7t=[  0.572  0.394  0.034] 

It  should  be  obvious  that  the  matrices  trained  directly  to  the  conceal  state(O).  The 
highest  probabilities  are  in  the  recurrent  conceal  state  and  transitions  to  the 
conceal  state  for  the  A  matrix.  The  B  matrix  obviously  favors  north,  n  favors  the 
conceal  state  initially. 


"0.572  0.613  0.497" 

y  =  0.394  0.365  0.464  The  expected  states  all  fall  to  the  conceal  state. 

0.034  0.022  0.039 


bestset  = 


The  highest  probabilities  favor  north  and  the  conceal  state.  This  is  obviously  so  for  the 
first  column. 
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Matrices  for  trivial  sequence  3,3,0 


■  0.488 

0.454 

0.059' 

■  0.463 

0.019 

0.019 

0.498  ■ 

A  = 

0.226 

0.577 

0.197 

B  = 

0.421 

0.019 

0.019 

0.541 

0.062 

0.46 

0.478 

0.23 

0.019 

0.019 

0.731 

71  =[0.022  0.33  0.649] 


It  should  be  obvious  that  the  matrices  trained  directly  to  the  speed  state{0).  The 
highest  probabilities  are  in  the  recurrent  speed  state  for  the  A  matrix.  The  B  matrix 
obviously  favors  the  observed  data,  ti  favors  the  speed  state  initially. 


‘  0.022 

0.029 

0.047' 

y  = 

0.33 

0.404 

0.571 

0.649 

0.567 

0.383 

The  expected  states  all  fall  to  the  speed  state  for  the 
first  two  iterations.  After  that,  some  indecision  is 
introduced  and  then  the  state  is  expected  to  go 
ambivalent.  The  results  make  sense  given  the  input 
data. 


bestset  = 


The  highest  probabilities  favor  west  and  the  speed  state  because  west  represents 
lowest  areas  of  the  field.  This  is  obviously  so  for  the  first  column  probability. 
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Matrices  for  trivial  sequence  0,1,2 


. !  0.561 

0.417 

0.022  ■ 

A=  0.239 

0.645 

0.116 

i  0.064 

0.572 

0.363 

It  =[0.502 

0.446 

0.052  ] 

■  0.414 

0.355 

0212 

0.02 

B  = 

0.279 

0.317 

0.385 

0.02 

0.209 

0.253 

0.518 

0.02 

It  should  be  obvious  that  the  matrices  trained  to  no  specific  state.  The  highest 
probabilities  are  in  the  recurrent  states  for  the  A  matrix.  The  B  matrix  obviously 
favors  the  observed  data.  «  doesn't  favor  any  particular  state  either. 


0.502 

0.43 

0.257  ■ 

0.446 

0.507 

0.615 

0.052 

0.063 

0.128 

The  expected  states  favor  conceal  and  ambivalent 
because  the  circle  it  made  headed  towards  the  conceal 
initially.  We  can't  really  pin  down  a  state  here  because 
there  was  no  definite  pattern.  The  object  would  stay  in 
the  ambivalent  state. 


bestset  = 


There  are  really  no  high  probabilities  here.  The  object  could  move  one  of  several 
ways.  This  is  more  evidence  that  the  object  is  ambivalent  as  to  were  he  wants  to  go. 


APPENDIX  J:  MODEL  ORDER 

First  we  will  lcx}k  at  the  mcxlels  with  4  observation  symbols  allowed. 


realO  := 


'1  2  2  2  2  l\ 
10  3  113 
l2  1  2  3  2  2] 


Real  ;=  5 


This  is  the  data  that  was  used  to  test  each  model.  The  LOGP 
routine  found  the  probability  of  each  sequence  from  left  to  right 
across  the  realO  matrix  Each  probability  was  converted  to  log 
form.  These  log  probabilities  are  the  output  of  the  model. 


real  :=0..  Real 

The  S  refers  to  the  number  of  states,  while  the  O  refers  to  the  number  of  observation  symbols  for 
the  model. 


- 1.174  ■ 

-.806  ■ 

- 1.297 

-.997 

-1.409 

- 1.438 

State3_Obs4  :  = 

State4_Obs4  := 

- 1.472 

- 1.585 

- 1.556 

-1.851 

-2.647 

-2.67 

Real 

Score_S3_04  =  State3_Obs4. 
i  =  0 

Score  S3  04  =-9.555 


Real 

Score_S4_04  =  State4_Obs4. 
i  =  0 


Score  S4  04  =-9.347 


-.873 


State5  Obs4  = 


1.085 

1.379 

1.413 

1.726 


-2.693 


Real 

Score_S5_04  =  State5_Obs4. 
i  =  0 


As  we  mentioned  before,  the  lowest  score 
will  show  the  most  appropriate  model  given 
the  realO  data. 


Score  S5  04  =-9.169 
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Score^  :=  Score_S3_04 

Scote^  :=  Score_S4_04 

Score^  :=  Score_S5_04  State  :=  3..  5 


Model  Score  versus  Allowed  State 


It  is  obvious  that  as  the  state  increases  the  model  becomes  more  appropriate.  Here  the  5  state 
model  would  be  best  because  it  has  the  lowest  score.  In  other  words,  the  5  state  model 
matrices  are  closest  to  being  ideal  so  we  should  use  that  model.  The  order  for  the  4 
observational  symbol  models  is  5. 
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Next  we  will  look  at  the  models  with  6  observation  symbols  allowed. 


realO  := 


1  2 
1  0 
1 


2  1 
1  1 
1  3 


2  l' 
2  3 
2  2, 


This  is  the  data  that  was  used  to  test  each  model.  The  LOGP 
routine  found  the  probability  of  each  sequence  from  left  to  right 
across  the  realO  matrix.  Each  probability  was  converted  to  log 
form.  These  log  probabilities  are  the  output  of  the  model. 


Real  =5 


real  0..  Real 

The  S  refers  to  the  number  of  states,  while  the  O  refers  to  the  number  of  observation  symbols  for 
the  model. 


- 1.847' 

- 1.583' 

-1.851 

- 1.589 

- 1.858 

- 1.734 

Statc3_Obs6  :  = 

Stat©4_Obs6  := 

- 1.874 

-1.93 

- 1.874 

- 1.997 

- 1.881 

-2.032 

Real 

Score_S3_06  .=  State3_Obs6. 
i  =  0 


Real 

Soore_S4_06  :=  State4_Obs6. 
i  =  0 


Score_S3_06  =-11.185 


Score_S4_06  =-10.865 


State5  Obs6  : 


- 1.527 
- 1.529 
-1.714 
- 1.903 
-2.04 


-2.086 


Real 

Score_S5_06  =  State5_Obs6. 
i  =  0 


As  we  mentioned  before,  the  lowest  score 
will  show  the  most  appropriate  model  given 
the  realO  data. 


Score  S5  06  =-10.799 
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Score  j  ;=Score_S3_06 

Score^  :=  Scoie_S4_06 

Scoie^  :=  Soon_SS_06  Stite  ;=  3..  S 


Model  Score  versus  Allowed  State 


It  is  obvious  that  as  the  state  irtcreases  the  model  becomes  more  appropriate.  Although  there 
is  a  difference  from  the  6  symbol  versus  the  4  symbol  model.  The  4  state  model  here  would 
probably  serve  just  as  well  as  the  5  state  model  because  there  is  such  a  small  change. 
Strictly  speaking,  however,  the  5  state  model  would  sGII  be  best 


